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(ABSTRACT) 


This  work  is  concerned  with  the  development  of  computational  methods  for  approximating 
sensitivities  of  solutions  to  boundary  value  problems.  We  focus  on  the  continuous  sensitiv¬ 
ity  equation  method  and  investigate  the  application  of  adaptive  meshing  and  smoothing 
projection  techniques  to  enhance  the  basic  scheme.  The  fundamental  ideas  are  first  devel¬ 
oped  for  a  one  dimensional  problem  and  then  extended  to  2-D  flow  problems  governed  by 
the  incompressible  Navier-Stokes  equations.  Numerical  experiments  are  conducted  to  test 
the  algorithms  and  to  investigate  the  benefits  of  adaptivity  and  smoothing. 
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Chapter  1 
Introduction 


During  the  past  decade,  considerable  effort  has  been  devoted  to  the  development  of  com¬ 
putational  tools  for  the  analysis,  design,  control  and  optimization  of  complex  physical 
systems.  Sensitivity  analysis  plays  an  important  role  in  all  aspects  of  this  effort.  Many 
simulation  tools  can  be  enhanced  by  using  local  information  provided  by  state  sensitivities 
to  obtain  near-by’  solutions.  For  example  Godfrey  ([20])  used  flow  sensitivities  to  reduce 
computational  time  in  simulations  of  multi-species,  chemically  reacting  flows.  Sensitivity 
methods  have  long  been  used  to  design  and  optimize  structures  (see  [23],  [18] ).  Moreover, 
sensitivity  reduction  is  often  employed  in  the  design  robust  feedback  controllers  (see  [15]’ 
[26],  [16],  [27],  [28]).  In  recent  years,  applications  of  gradient  based  optimization  algorithms 
to  multi-disciplinary  design  optimization  (MDO)  problems  have  lead  to  renewed  interest 
in  state  sensitivities  as  one  method  to  compute  cost  function  gradients  (see  [4],  [7]). 

In  order  to  make  effective  use  of  sensitivities  in  analysis,  design,  and  control,  accurate  and 
efficient  computational  algorithms  are  essential.  In  addition,  the  natural  trade-off  between 
accuracy  and  efficiency  can  often  be  exploited  to  produce  the  best  overall  computational 
tool  for  a  specific  application.  For  example,  when  used  to  enhance  a  simulation  tool,  accu¬ 
racy  may  be  more  important  than  speed.  However,  when  used  for  gradient  computations 
in  an  optimization  loop,  speed  may  be  more  crucial  than  high  accuracy. 

Two  early  numerical  techniques  employed  for  sensitivity  calculations  were  the  straight¬ 
forward  finite  difference  calculations  and  the  “discretize-then-differentiate”  approach.  The 
“discretize-then-differentiate”  scheme  approximates  sensitivities  by  first  employing  some 
discretization  scheme  to  approximate  the  solution  to  a  PDE  and  then  implicitly  differenti¬ 
ating  this  result  to  to  obtain  a  sensitivity  approximation  scheme.  Recent  advances  in  this 
area  include  the  development  of  automatic  differentiation  packages  (see  [32],  [21]),  like  AD- 
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IFOR,  which  given  source  code  for  state  calculations,  generates  source  code  for  sensitivity 
calculations. 

There  are  a  number  of  significant  disadvantages  to  each  of  the  above  approaches.  One 
disadvantage  of  finite  difference  techniques  is  that  they  are  computationally  intensive  since 
they  require  two  flow  solves  at  different  parameter  values  to  form  a  difference  quotient.  In 
large  aerospace  flow  problems,  this  is  sometimes  an  extremely  severe  requirement  because 
mesh  generation  itself  can  take  weeks  to  months.  In  addition,  it  is  not  generally  known 
a-priori  what  is  a  sufficiently  small  step  size  to  ensure  accurate  approximations,  making  it 
difficult  to  incorporate  these  ideas  into  optimizations  schemes.  A  procedure  for  estimat¬ 
ing  an  optimal  step-size  is  described  in  [30].  The  optimal  step-size  attempts  to  balance 
truncation  error  and  round-off  error.  Unfortunately,  since  the  discretization  is  generally 
non-uniform,  the  truncation  error  varies  from  point  to  point  so  that  the  results  must  be 
used  in  combination  with  some  other  decision  algorithm. 

A  disadvantage  of  the  discretize-then-differentiate”  technique  is  that  in  cases  where  the 
mesh  is  parameter  dependent,  as  is  the  case  in  shape  optimization  problems,  then  differ¬ 
entiation  of  the  discrete  PDE  leads  to  mesh  sensitivities  on  the  right  hand  side.  Although 
there  has  been  much  work  done  in  recent  years  to  get  a  handle  on  these  quantities  (see 
[38]),  calculating  mesh  derivatives  is  still  not  well  understood,  particularly  in  cases  where 
the  meshes  are  prescribed  adaptively. 

In  recent  years  several  new  approaches  have  been  developed  in  an  attempt  to  alleviate 
some  of  the  disadvantages  of  the  two  approaches  mentioned  above.  The  idea  is  to  em¬ 
ploy  a  “differentiate-then-discretize”  scheme  or  the  so-called  continuous  sensitivity  equa¬ 
tion  method  (SEM).  The  SEM  consists  of  implicitly  differentiating  the  PDE  to  obtain  a 
sensitivity  equation  (SE).  Then  both  the  PDE  and  the  SE  are  discretized  to  obtain  finite 
dimensional  equations  for  numerical  approximations  to  the  PDE  and  the  SE. 

This  technique  reduces  some  of  the  problems  mentioned  above.  The  SEM  eliminates  the 
need  for  a  second  flow  solve  saving  costly  computer  time.  Also,  since  the  differentiation  is 
done  before  the  discretization,  it  eliminates  the  need  to  calculate  “mesh”  derivatives,  an 
exceedingly  difficult  task  especially  in  the  case  of  complicated  2-D  and  3-D  flow  problems. 
The  method  has  another  advantage.  The  sensitivity  equations  are  always  linear  and,  in 
principal,  it  is  simple  to  modify  existing  flow  solvers  so  that  one  additional  Newton  step 
at  the  end  of  a  flow  solve  is  sufficient  to  approximate  sensitivities.  This  feature  makes  the 
SEM  an  easy  approach  to  implement. 


The  SEM  produces  an  approximation  for  the  solution  to  the  continuous  sensitivity  equation. 
In  some  cases,  especially  in  the  case  of  optimization,  what  one  really  needs  is  the  sensitivity 
of  the  discrete  solution.  These  solutions  are  not  necessarily  equivalent.  Because  of  this 
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discrepancy,  it  is  not  obvious  that  the  SEM  would  lead  to  convergence  of  the  optimization 
algorithm.  However,  recent  theoretical  and  numerical  studies  have  shown  that  this  issue 
can  often  be  addressed. 

In  [11],  Burkardt  applied  this  sensitivity  approximation  scheme  to  a  shape  optimization 
problem  involving  incompressible  fluids.  He  found  that  in  many  cases  the  SEM  gave  suf¬ 
ficiently  good  approximations  to  allow  the  optimization  algorithm  to  converge.  However, 
Burkardt  also  provided  examples  where  the  differences  between  the  discrete  sensitivities 
and  the  numerical  solution  of  the  continuous  sensitivity  equation  were  large  enough  to  de¬ 
stroy  the  convergence  of  the  optimizer.  This  was  true  even  in  cases  where  the  finite  element 
method  had  done  a  good  job  approximating  the  solution  of  the  continuous  problem.  He 
also  found  some  problems  at  higher  Reynolds  number  flows.  One  of  the  issues  which  may 
have  been  a  factor  is  the  accuracy  of  the  flow  derivative  information  used  in  the  sensitivity 
calculation.  We  will  show  in  this  paper  that  improved  derivative  information  can  dras¬ 
tically  improve  the  sensitivity  approximations  and,  in  turn,  have  dramatic  effects  on  the 
performance  of  an  optimization  scheme. 

In  [3],  Borggaard  provided  a  partial  analysis  of  this  problem.  He  proved  that  if  the  discrete 
approximation  to  the  continuous  sensitivity  equation  was  “close  enough”  to  the  sensitiv¬ 
ity  of  the  discrete  solution,  then  an  appropriately  chosen  optimization  algorithm  would 
converge.  He  also  showed  that  using  the  SEM,  rather  than  the  finite  difference  approach, 
could  reduce  CPU  times  by  50  percent  or  more.  This  is  an  extremely  important  savings 
for  high-cost  flow  problems. 

The  influence  of  flow  discontinuities  on  sensitivity  approximations  was  investigated  by 
Appel  in  [2].  Appel  compared  all  three  sensitivity  approximation  techniques  on  Euler  flows 
as  well  as  the  1-D  Riemann  problem,  a  problem  for  which  an  exact  solution  is  known.  He 
found  that  the  shocks  in  the  solution  led  to  errors  in  the  sensitivity  approximations  for  all 
of  the  techniques. 

We  are  interested  in  applications  that  lend  themselves  to  physics  based  modeling.  Specif¬ 
ically,  we  shall  concentrate  on  problems  where  the  governing  equations  are  differential 
equations.  We  use  a  simple  boundary  value  problem  in  one  spatial  dimension  (1-D)  to  in¬ 
troduce  the  ideas  and  illustrate  the  various  methods.  This  model  problem  is  described  by 
an  ordinary  differential  equation.  Once  the  basic  ideas  and  techniques  have  been  tested  on 
this  1-D  problem,  we  turn  to  more  substantial  problems  from  fluid  dynamics.  In  particular, 
we  consider  2-D  fluid  flows  governed  by  the  Navier-Stokes  equations. 


Chapter  2 

A  1-D  Model  Problem 


In  this  chapter,  we  consider  a  simple  1-D  non-linear  boundary  value  problem.  We  formulate 
an  optimization  problem  and  provide  a  detailed  presentation  of  the  SEM.  We  use  finite 
elements  to  construct  approximations  and  compare  the  numerical  results  to  exact  solutions. 
The  goal  of  this  chapter  is  to  use  the  model  problem  to  describe  the  basic  ideas,  identify 
the  fundamental  issues,  and  set  the  stage  for  more  complex  problems  to  come. 


2.1  Model  Problem 


We  concentrate  on  the  boundary  value  problem  defined  by  the  non-linear  differential  equa¬ 
tion 


cP  ,  ,  1 

^w(*>  +  8 


*«<*) 


=  0  for  0  <  x  <  q, 


with  boundary  conditions 

w(0)  -  0  and  w(q)  =  4. 

For  each  q  >  1,  the  exact  solution  to  this  boundary  value  problem  is  given  by 


w(x)  =  w(x,q)  =  i\jx+  — - 2(g-  1). 


(2.1) 


(2.2) 


(2.3) 


Let  0  <  X\  <  £2  <  ...  <  xp  <  1  be  fixed  locations  and  assume  that  Wj  is  data  representing 
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values  of  w(x)  at  Xj.  Consider  the  inverse  design  problem:  Find  q*  >  1  such  that 

Jp(q*)  <  Jp(q )  =  i  q)  -  Wj))2,  (2.4) 

3= 1 

where  w(x,g)  is  the  solution  of  (2.1)  -  (2.2)  and  the  integer  p  represents  the  number  of 
data  points.  In  gradient-based  optimization  one  needs  the  derivative 

P  Q 

VqJM  =  ?)  -  Wj)—w(xj,  q).  (2.5) 

One  approach  to  the  evaluation  of  this  gradient  at  q  is  to  “compute”  the  state,  w(xj,  q), 
the  sensitivity,  -§^w(xj,q),  and  form  the  computation  (2.5).  This  involves  first  solving  (2.1) 
-  (2.2)  for  w(x,q )  and  then  computing  the  sensitivity  -^w(x,q). 


2.2  The  Sensitivity  Equation 


One  benefit  of  using  the  model  problem  is  that  we  can  calculate  the  sensitivity  -§-w(x,  q) 
by  direct  differentiation  of  (2.3).  In  particular,  9 


d  f  N  9-1 

-W(x,g)  =  - - 


-  2. 


'x  + 


(<7~1)2 


(2.6) 


On  the  other  hand,  we  can  implicitly  “differentiate”  the  boundary  value  problem  (2.1)  - 
(2.2)  and  obtain  a  boundary  value  problem  for  the  sensitivity  -§^w(x,q)  =  s(x,q )  =  s(x). 
It  follows  that  s(x)  satisfies  the  linear  differential  equation 


d2  3 

+  8 


d  ,  , 

SW(X) 


d_ 

dx 


with  boundary  conditions 


s(x)  =  0, 

(2.7) 

9  t  \ 

~ Txw(q) ■ 

(2.8) 

It  is  important  to  be  cautious  when  using  (2.2)  to  derive  the  boundary  conditions  for  the 
sensitivity  equation.  The  first  boundary  condition, 


Dawn  L.  Stewart 


Chapter  2.  A  1-D  Model  Problem 


6 


is  rather  obvious.  However,  deriving  the  second  boundary  condition  in  (2.8)  can  be  tricky 
(especially  in  2-D  and  3-D  domain  optimization  problems).  Using  the  boundary  condition 
(2.2),  the  chain  rule  leads  to  the  correct  boundary  condition 


0=  ~w 
Dq 


(g,q)  =  {JU(x,9)  +  -j^w(x,q) } 


Since  this  is  a  model  for  a  more  difficult  applied  problem  for  which  we  would  not  know  the 


state,  w,  and  the  sensitivity,  s,  exactly,  we  are  interested  in  constructing  good  methods  for 
obtaining  numerical  approximations  of  these  quantities. 


2.3  Domain  Mapping 


For  computational  purposes,  a  series  of  transformations  are  used  to  first  map  the  prob¬ 
lem  to  a  fixed  “computational  domain”  and  second  to  transform  the  non-homogeneous 
boundary  conditions  to  homogeneous  Dirichlet  conditions.  This  is  standard  in  many  CFD 
(computational  fluid  dynamics)  algorithms.  What  is  not  typical  is  that  we  apply  a  similar 
transformation  to  the  sensitivity  equation.  Let  T  :  [0,  <j\  ->  [0, 1]  be  defined  by 

T(x,q)  =  ^=C  (2.9) 

To  avoid  confusion,  we  use  x  for  the  independent  variable  on  [0,  q]  and  £  for  the  independent 
variable  on  [0,1].  Note  that  T  =  T(x,q )  depends  explicitly  on  the  “shape  parameter”  q. 
For  each  q  >  1,  the  inverse  transformation 

\T(x,  ?)]  1  =  M (£,  q) :  [0, 1]  — >•  [0,  g] 

is  defined  by 


M(£,  q)  =q£  =  x. 


(2.10) 


Let  w(g)  =  w(C  q)  =  w(M(£,q),q)  and  define  z(£)  by 

z(Z)  =  q)  =  w(0  -  4f.  (2.11) 

Applying  this  transformation,  we  obtain  the  Dirichlet  problem  defined  on  the  computational 
domain  [0, 1]  by  the  differential  equation 


~iez(C 1  +  % 


d_ 


*(0+4 


=  0 


(2.12) 
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with  Dirichlet  boundary  conditions 

z(0)  =  0  and  z(  1)  =  0.  (2.13) 


In  a  similar  manner,  we  transform  the  sensitivity  equation  (2.7)  -  (2.8)  to  the  computational 
domain  [0, 1].  In  particular,  we  let  s(£)  =  s(£,  q)  =  s(M(£,q),q)  and  define  p(£,  q)  by 

— u)(l) 

p(0  =  p(£>  q)  =  s{0  +  — £.  (2.14) 

It  follows  that  the  transformed  sensitivity  p(£)  satisfies  the  equation 


d2  3 

—p(0  +  — 


d e: 


8,  l!*®*4’ 


d  j|z(l)  +  4' 


A 


p(0  - 


with  Dirichlet  boundary  conditions 

p(o)  =  0 


and  p(l)  =  0. 


(2.15) 


(2.16) 


In  practice,  one  must  use  some  numerical  scheme  to  solve  the  boundary  value  problem  (2.1) 
-  (2.2),  and  the  computation  of  -§^w(x,  q)  must  be  accomplished  by  using  this  approximate 
solution.  We  approach  this  problem  by  solving  the  corresponding  sensitivity  equation. 
As  shown  below,  there  are  many  “natural”  numerical  schemes  that  one  can  employ  in  this 
approach.  Although  we  discuss  several  schemes,  we  will  concentrate  on  a  projection  method 
approach  in  later  chapters.  The  basic  idea  can  be  extended  to  complex  aerodynamic  flow 
problems.  However,  many  theoretical  and  technical  issues  are  not  yet  settled. 

Comment:  It  is  important  to  note  that  the  construction  of  the  transformed  state  equation 
(2.12)  -  (2.13)  and  the  transformed  sensitivity  equation  (2.15)  -  (2.16)  requires  the  derivative 
(in  space)  of  the  transformation  In  particular,  one  needs  ^  M(£,q)  or  else  a 

numerical  approximation  of  %M($,q).  This  issue  is  addressed  in  many  CFD  codes  and 
there  are  good  methods  for  dealing  with  this  problem  (e.g.  see  [38]).  However,  there  is  no 
need  to  compute  the  derivative  -fcMfcq)  with  this  approach.  On  the  other  hand,  if  one 
transforms  the  state  equations  and  then  derives  the  sensitivity  equation  for  the  transformed 
state  equation,  then  the  chain  rule  requires  the  calculation  of  j-M(£,q).  In  particular,  if 
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d{£>q)  then 

d{£,q)  =  q) 

=  -^w(M(^q),q) 

$  $  d 

=  —w(M  (£,  q),q)  +  q ),  g)^M(£,  ?) 

^  o 

=  s&q)  +  -^w(M(Z,q),q)—M(t,q). 

Consequently,  the  difference  between  the  transformed  sensitivity  s(£)  and  the  sensitivity  of 
the  transformed  state  d(f)  is  given  by  the  difference 

q)  -  ■*(£,  q)  =  -^w(M(C  q),  q)^M(^  q).  (2.17) 

The  right  hand  side  of  (2.17)  is  the  continuous  version  of  the  “mesh”  gradient  and  the 
source  of  considerable  computational  complexity.  Consequently,  one  advantage  of  mapping 
both  the  state  and  the  sensitivity  equation  is  that  the  computation  of  this  gradient  can  be 
eliminated. 

Finally,  we  note  that  once  z(£,q)  and  p(£,q)  are  computed  on  [0, 1],  the  state  w(x,q )  and 
sensitivity  s(x,  q)  can  be  recovered  on  [0, g]  by 

w(x,  q)  =  w(T(x,  q),  q)  =  z(T(x,  q),q)  +  4 T(x,  q),  (2.18) 

and 

s(a;,  q)  =  s(T(x,  q),  q)  =  p(T(x,  q),  q)  -  +  -  T(x,  q),  (2.19) 

respectively.  When  applying  the  inverse  transforms  to  the  numerical  solutions,  there  is  the 
possibility  that  numerical  errors  can  be  induced.  In  particular,  the  map  T(x,  q)  is  often 
constructed  numerically  for  practical  CFD  problems.  Moreover,  in  equation  (2.19)  the 
presence  of  the  derivative  ^z(  1,  q)  at  the  boundary  can  introduce  additional  errors.  These 
are  practical  issues  that  are  important  to  address  in  more  complex  problems.  However, 
observe  once  again  that  it  is  not  necessary  to  compute  a  mesh  gradient  to  transform  back 
to  the  physical  domain. 
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2.4  Computational  Algorithms 


We  turn  now  to  the  issue  of  approximating  the  coupled  system  defined  by  the  state  equation 


d2  1 

sr(f)+i; 


1  3 


k 


z(£)  +  4 


=  0 


(2.20) 


and  sensitivity  equation 


d 2  _  3 

—p(0  +  ~ 


<%■ 


8  q 


1  2 


^(0  +  4 


+4^ 


=  0 


(2.21) 


with  boundary  conditions 


2(0)  =  2(1)  =  0, 

and 


(2.22) 


p(0)  =  p(l)  =  0,  (2.23) 

respectively.  It  is  important  to  observe  that  (2.20)-(2.21)  is  weakly  coupled  in  the  sense 
that  (2.21)  does  not  feed  back  into  (2.20).  We  take  advantage  of  this  structure  to  develop 
a  family  of  numerical  schemes  for  computing  the  sensitivity  p(f). 

Note  that  for  each  q  >  1,  z(£,q)  and  p(f ,  q)  exist  and  belong  to  J3j(0,  l)f)^2(0, 1).  Let 
V  =  i?o(0, 1)  x  Ifj (0,1)  and  observe  that  for  each  (<£(-),  ^(-))  €  V  the  solution  pair 
(z(-),p(-))  satisfies  the  weak  system 


(z  O)’0  (•)}  +  (•)  +  4  ,<K-)^=  0 


-  (p  (■),*' (■))  +  ~  ([*'(•)  +4]2  (p(.)  - 
where  (•,  •)  denotes  the  L2  inner-product. 


+  4 


=  0 

(2.24) 

(2.25) 

2.5  A  Finite  Element  Scheme 


Although  there  are  several  possible  choices  for  finite  element  spaces,  we  shall  limit  our 
discussion  to  the  simplest  (convergent)  scheme.  We  note  that  in  more  complex  problems 
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one  must  choose  these  spaces  with  care  to  ensure  the  algorithm  satisfies  the  appropriate 
convergence  criteria  (inf-sup  conditions,  etc.).  Let  h  =  and  denote  by  Vh  the 

product  space 

F‘  =  ftl)xftl)CV,  (2.26) 

where  «Sf  ( 0, 1)  denotes  the  subspace  of  i70x(0, 1)  consisting  of  continuous  piecewise  linear 
functions  with  nodes  at  &  =  i  =  1, 2, . . .  ,  K.  Let  /if  (•)  and  hf  (•)  denote  the  standard 
“hat”  functions  with  nodes  &  and  respectively.  Also  let 

AO  =5^7if(£)  (2.27) 

t=l 


and 


M 


pM(« = (?) 

i= i 


(2.28) 


be  Galerkin  approximations  of  the  pair  ( z(-),p( •))  in  V\  Consider  the  Galerkin  approxi¬ 
mations 


<•))+£ 


r"(-)+4 


(2.29) 


(2.30) 


Note  that  pM(-)  depends  on  **(•)  and  its  spatial  derivative  -%zN(-)  on  [0, 1],  To  emphasize 

the  dependence  we  let  pN,M (•)  denote  the  solution  of  (2.30),  given  that  zN (•)  obtained  from 
(2.29)  is  used  in  (2.30),  and  let 


vh(^q)  =  (zN(Z,q),pN’M(Z,q))  (2.31) 

denote  the  solution  pair.  At  this  point,  there  are  two  important  observations  that  play  a 
key  role  in  the  construction  of  accurate  numerical  sensitivities. 
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•  The  freedom  to  choose  separate  finite  element  spaces  for  the  state  z(£,q)  and  the 
sensitivity  pfa  q)  allows  for  the  development  of  schemes  that  simultaneously  converge 
to  the  state  and  sensitivity.  In  addition,  both  h-refinement  (mesh  refinement)  and 
p-refinement  (the  selection  of  higher  order  elements)  can  be  combined  to  construct 
numerical  solutions  of  vh(£,q)  =  (zN(£,  q),pN’M({,q))  so  that  the  error  in  pN’M(£,q) 
is  sufficiently  small  to  ensure  convergence  of  optimal  design  algorithms  based  on  the 
SEM  (see  [3],  [7],  [5]). 

•  The  solution  pN’M(£)  depends  not  only  on  zN(£)  but  also  its  derivative  -^zN(t). 
Moreover,  since  zN(£)  is  piecewise  linear,  ± zN( f)  is  a  piecewise  constant  (discontin¬ 
uous)  function.  However,  the  actual  transformed  sensitivity  p({)  is  smooth,  and  one 
might  expect  to  lose  at  least  one  order  of  accuracy  in  pN’M({).  In  fact,  things  can  be 
much  worse  unless  special  care  is  exercised. 


There  are  two  obvious  “fixes”  to  address  these  issues.  One  could  use  higher  order  splines 
for  the  sensitivity  variable  p(£).  However,  this  method  will  be  more  expensive,  and  it  is 
not  reasonable  to  expect  great  improvements  unless  higher  order  schemes  are  also  used 
for  the  state  equation.  The  other  obvious  fix  is  to  use  mesh  refinement  in  M  (assuming 
accuracy  in  N).  There  is  a  third  approach  that  makes  use  of  “smoothing  projections.”  The 
idea  is  similar  to  the  method  used  to  obtain  a-posteriori  error  estimators  for  adaptive  mesh 
generation  (see  [8],  [24],  [41],  [42]).  This  approach  will  be  outlined  in  Chapter  4  and  applied 
to  the  model  problem  and  to  a  2-D  fluid  flow  problem. 


2.6  Numerical  Results 


We  use  the  FEM  scheme  outlined  above  to  construct  approximations  to  the  state  and  the 
state  sensitivities.  We  also  evaluate  the  use  of  the  state  sensitivities  in  an  optimization 
algorithm  to  solve  the  inverse  design  problem  outlined  in  Section  2.1. 


2.6.1  Convergence  of  Solutions  for  the  Boundary  Value  Problem 
and  Sensitivity  Equation 

In  this  section,  we  compare  the  finite  element  approximations  of  the  solutions  to  the  state 
and  sensitivity  equations  with  their  exact  solutions.  First,  we  note  that  the  finite  element 
scheme  converges  to  the  exact  solution  of  the  nonlinear  problem  (for  each  q  >  1).  Figure 
(2.1)  shows  the  finite  element  approximations  to  the  solution  of  the  boundary  value  problem 


Dawn  L.  Stewart 


Chapter  2.  A  1-D  Model  Problem 


12 


Figure  2.1:  Numerical  Approximations  to  the  Solution  of  the  Boundary  Value  Problem  at 
q  =  2  and  q  =  1.2 


at  two  parameter  values:  q  =  2  and  q  =  1.2.  Notice  that  at  q  =  2,  the  N  =  4  finite  element 
model  provides  an  excellent  match  to  the  exact  solution.  However,  when  q  =  1.2  one 
sees  that  a  finer  mesh  (N  =  8)  is  required  to  obtain  the  same  order  of  accuracy.  This 
convergence  is  expected  because  the  gradient  of  the  solution  becomes  singular  as  q  -*  1+ 
and  hence  the  problem  becomes  stiff  in  this  parameter  region.  This  is  also  the  case  for  the 
sensitivity  equation. 

Consider  the  corresponding  finite  element  approximations  of  the  sensitivity  equation. 
Recall  that  N  and  M  define  the  meshes  for  the  state  and  sensitivity  equations,  respectively. 
Figure  (2.2)  shows  the  finite  element  approximations  for  the  sensitivity  with  q  =  2  and 
N  =  M  ranging  from  2  to  16.  Observe  in  Figure  (2.3)  that  although  the  finite  element 
scheme  produces  excellent  solutions  to  the  state  equation  when  N  =  4,  the  error  in  the 
corresponding  sensitivity  does  not  diminish  until  N  =  M  =  16. 

As  mentioned  before,  one  obvious  “fix”  is  to  use  mesh  refinement  in  M.  Figure  (2.4) 
shows  the  results  for  this  technique  when  N  =  2  and  M  ranges  from  2  to  16.  Note  that 
improvements  in  the  accuracy  of  the  sensitivity  approximation  are  limited  by  the  accuracy 
of  the  state  approximation.  When  we  increased  N  to  4,  the  sensitivity  errors  were  decreased 
(see  Figure  (2.5)). 

The  stiffness  of  the  problem  near  q  =  1  increases  the  difficulty  of  getting  good  sensitivity 
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Figure  2.2:  Numerical  Approximations  to  the  Solution  of  the  Sensitivity  Equation  at  q  =  2 
using  piecewise  constant  (PWC)  derivatives 


- j - r - , 

T-  , 

Lrror  for  flow 

.  Error  for  sensitivity 
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Figure  2.3:  L ?  Error  of  the  Solution  and  Sensitivity  Approximations  at  q  =  2 
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Figure  2.4:  Numerical  Approximations  to  the  Solution  of  the  Sensitivity  Equation  at  q  =  2 
using  PWC  derivatives  with  N  =  2  and  mesh  refinement  in  M 

approximations.  Figures  (2.6)  and  (2.7)  show  that  the  sensitivity  approximations  become 
unreliable  as  q  — »•  1+  even  though  the  state  approximations  are  still  fairly  good.  Figure 

(2.8)  displays  the  difference  between  the  L 2  error  in  the  state  approximation  and  the  L2 
error  in  the  sensitivity  approximation.  As  N,  M  increases  we  obtain  convergence  of  the 
scheme,  but  the  approximations  for  smaller  N ,  M  contain  large  errors  and  the  convergence 
of  the  finite  element  approximation  to  the  analytical  solution  is  not  at  all  monotone.  Figure 

(2.9)  is  a  graph  of  the  L2  error  of  the  sensitivity  approximations  for  various  values  of  q. 
We  shall  observe  similar  behavior  in  the  next  section,  where  we  approximate  solutions  to 
2-D  flow  problems. 
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q  =2  a  =? 


Figure  2.5:  Numerical  Approximations  to  the  Solution  of  the  Sensitivity  Equation  at  q  =  2 
using  PWC  derivatives  with  N  =  4  and  mesh  refinement  in  M 


X 


Figure  2.6:  Numerical  Approximations  to  the  Solution  of  the  Sensitivity  Equation  at  q  =  1.4 
using  piecewise  constant  (PWC)  derivatives 
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Figure  2.7:  Numerical  Approximations  to  the  Solution  of  the  Sensitivity  Equation  at  q  =  1.2 
using  piecewise  constant  (PWC)  derivatives 


Figure  2.8:  L 2  Error  of  the  Solution  and  Sensitivity  Approximations  at  q  =  1.2 
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Figure  2.9:  L2  Error  of  Sensitivity  Approximations  using  Piecewise  Constant  (PWC) 
Derivatives 
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2.6.2  Optimization  Results 


We  now  return  to  the  inverse  design  problem  presented  in  §2.1  and  evaluate  the  conver¬ 
gence  properties  of  an  optimization  scheme  using  the  sensitivity  equation  method.  After 
discretization,  the  infinite  dimensional  inverse  problem  (2.4)  becomes-  Find  a*  >  1  such 
that 


Jp  (?*)  ^  Jp  (?)  =  |  ^2(wN(xj:  q )  -  wj)2, 


(2.32) 


3= 1 


where  wN(x,q)  is  obtained  using  (2.18).  Notice  that  the  gradient  has  the  form 

= !>"(%,?)  - 

3= 1  " 


(2.33) 


The  sensitivity  equation  method  applied  to  the  optimization  problem  replaces  (£j,q), 
the  sensitivity  of  the  discrete  solution,  with  an  approximation  to  s(i,-,  o),  sN’M(xi  q)  from 


The  standard  Gauss-Newton  algorithm  is  used  to  approximate  q*  .  The  algorithm  solves  a 
least  squares  problem  at  each  iteration  and  proceeds  as  described  in  Table  2.1: 


Given  q0  and  Wi 

Set  iteration  counter:  i  =  0 

Compute  wf  =  wN(qi) 

Compute  sf’M  =  sN>M(qi) 

While  ( i  <  max  iterations  AND  ||VgJAr(gi)||  >  tolerance  )  Do 
Calculate  step  Sq 
Set  ql+i  =  qt  +  Sq 
Compute  wf+1  =  wN(qi+1 ) 

Compute  4’f  =  sN'M(qt+1) 

Increment  Counter 
EndWhile 


Table  2.1:  Gauss-Newton  Algorithm 


The  data  to  be  matched,  denoted  by  tbj  in  (2.4),  is  indicated  by  pluses  in  Figures  (2.10) 
-  (2.11).  This  data  set  was  generated  by  randomly  perturbing  the  value  of  w(xi)  in  (2.3) 
using  q  =  2  and  q  =  1.4  with  p  =  4,16  data  points.  Table  2.2  shows  the  numerical  values 
of  the  data  for  comparison  purposes. 
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P— 4 

q= 

2 

q=1.4  1 

Xj 

w{Xj) 

Wj 

perturbation 

w(Xj) 

Wj 

perturbation 

0.1250 

0.4495 

0.4516 

0.0021 

0.8248 

0.7431 

-0.0817 

0.2500 

0.8284 

0.7722 

-0.0563 

1.3541 

1.3875 

0.0335 

0.5000 

1.4641 

1.4080 

-0.0561 

2.1394 

1.9848 

-0.1546 

0.7500 

2.0000 

1.8338 

-0.1662 

2.7553 

2.5758 

-0.1795 

II 

Oi 

q= 

2 

q=1.4 

Xj 

w(£j ) 

Wj 

perturbation 

w(Xj) 

Wj 

perturbation 

0.0312 

0.1213 

0.1126 

-0.0088 

0.2677 

0.2690 

0.0013 

0.0625 

0.2361 

0.2207 

-0.0154 

0.4806 

0.4480 

-0.0326 

0.0938 

0.3452 

0.3670 

0.0218 

0.6629 

0.6375 

-0.0254 

0.1250 

0.4495 

0.4372 

-0.0123 

0.8248 

0.7563 

-0.0685 

0.1562 

0.5495 

0.5916 

0.0421 

0.9720 

0.9772 

0.0052 

0.1875 

0.6458 

0.6611 

0.0154 

1.1079 

1.0335 

-0.0744 

0.2188 

0.7386 

0.7651 

0.0265 

1.2347 

1.2820 

0.0473 

0.2500 

0.8284 

0.8687 

0.0403 

1.3541 

1.4106 

0.0565 

0.3125 

1.0000 

0.9091 

-0.0909 

1.5749 

1.4821 

-0.0928 

0.3750 

1.1623 

1.2674 

0.1051 

1.7768 

1.6265 

-0.1503 

0.4375 

1.3166 

1.2486 

-0.0680 

1.9641 

1.9460 

-0.0181 

0.5000 

1.4641 

1.6086 

0.1445 

2.1394 

2.0504 

-0.0890 

0.5625 

1.6056 

1.6483 

0.0427 

2.3048 

2.0950 

-0.2098 

0.6250 

1.7417 

1.9067 

0.1651 

2.4619 

2.3973 

-0.0646 

0.6875 

1.8730 

1.9668 

0.0938 

2.6117 

2.8040 

0.1922 

0.7500 

2.0000 

2.0554 

0.0554 

2.7553 

2.9296 

0.1743 

Table  2.2:  Matching  Data  for  the  Optimization 
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Figure  2.10:  Data  generated  at  q  =  2 


Figure  2.11:  Data  generated  at  q  =  1.4 
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The  exact  cost  functional,  J(q),  and  several  approximations  to  it,  JN(q)  are  plotted  in 
Figures  (2.12)-(2.13)  for  the  case  where  p  —  16  and  q*  ~  2  and  q*  ~  1.4,  respectively. 

Although  we  made  several  runs  with  various  data  sets,  we  present  the  results  of  two  runs. 

•  Case  1:  The  solution  was  calculated  for  q  =  2.0  and  the  noise  vector  in  Table  (2.2) 
was  added  to  obtain  data  for  optimization.  Here,  the  optimal  q  is  approximately 
q*  ~  2.  The  optimization  algorithm  was  started  at  qinit  =  1.2. 

•  Case  2.  In  this  case,  the  solution  was  calculated  for  q  —  1.4  and  the  noise  vector  in 
Table  (2.2)  was  again  added  to  obtain  data  for  optimization.  Here,  the  optimal  q  is 
approximately  q*  ~  1.4.  Here,  the  optimization  algorithm  was  started  at  qinit  =  2.0. 


The  scheme  was  considered  converged  when  the  norm  of  the  gradient  of  the  cost  functional 
was  less  than  10  .  Notice  that  the  simulations  were  performed  using  sensitivities  calculated 
using  the  natural  piecewise  constant  finite  element  gradient  approximations.  Tables  (2.3)  - 
(2.4)  show  the  results  of  these  simulations  as  N ,  M  ranged  from  2  to  128.  The  time  for  the 
runs  was  measured  in  seconds  and  the  runs  were  performed  on  a  Silicon  Graphics  Onyx2. 
Notice,  the  effect  of  the  bad  sensitivities  on  the  convergence  of  the  optimization  scheme  for 
Case  1  when  N  =  M .  As  expected,  larger  values  of  N,  M  were  required  for  convergence 
of  the  optimization  algorithm  for  Case  2. 


In  Chapter  4,  we  return  to  this  problem  and  use  smoothing  projections  to  enhance  sen¬ 
sitivity  computations  and  convergence.  We  turn  now  to  the  2-D  Navier-Stokes  equations 
and  discuss  two  specific  flow  problems. 
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N 


2 

2 

4 

4 

8 

8 

16 

16 

32 

32 

64 

64 

128 


M 


2 

5 

4 

9 

8 

17 

16 

33 

32 

65 

64 

129 

128 


CONV/DNC 

ITERATIONS 

TIME 

ii.ni 

q 

DNC 

20 

31.4 

2.4088 

1.21993 

CONV 

17 

3.1 

0.32610 

1.94862 

DNC 

20 

68.3 

2.6252 

1.15545 

CONV 

15 

6.6 

0.26469 

1.94795 

DNC 

20 

143.1 

3.1942 

1.08641 

CONV 

11 

7.4 

0.24332 

1.94328 

DNC 

20 

334.7 

3.0308 

1.16909 

CONV 

12 

16.0 

0.24235 

1.93902 

DNC 

20 

916.6 

3.1796 

1.18975 

CONV 

20 

42.1 

0.24278 

1.93771 

DNC 

20 

4323.1 

0.24276 

1.93477 

CONV 

12 

191.9 

0.242783 

1.93747 

CONV 

10 

252.2 

0.24285 

1.93584 

Table  2.3:  Optimization  Results  for  p  =  16,  q m  ~  2,  qinit  =  1.2,  PWC  Derivatives 


N 


2 

4 

4 

8 

8 

16 

16 

32 

32 

64 

64 

128 


M 


2 

5 

4 

9 

8 

17 

16 

33 

32 

65 

64 

129 

128 


CONV/DNC 


DNC 

DNC 

DNC 

DNC 

DNC 

DNC 

DNC 

CONV 

CONV 

CONV 

CONV 

CONV 

CONV 


ITERATIONS 


20 

20 

20 

20 

20 

20 

20 

18 

11 

16 

13 

17 

12 


TIME 


0.6 

32.5 

63.3 

68.4 
134.9 

61.3 

293.2 

18.8 

6.8 

13.2 

33.7 

58.0 

192.7 


II J 


W] 


0.86393 

0.72436 

0.54447 

0.52754 

0.60045 

0.50601 

0.78813 

0.41994 

0.41001 

0.41046 

0.40820 

0.40858 

0.40774 


1.83807 

1.33907 

1.48191 

1.4419 

1.38030 

1.47419 

1.33031 

1.43638 

1.41986 

1.42342 

1.41681 

1.42048 

1.41605 


Table  2.4:  Optimization  Results  for  p  =  16,  qm  ~  1.4,  qinit  =  2.0,  PWC  Derivatives 


Chapter  3 

2-D  Flow  Problems 


We  now  turn  to  problems  in  fluid  dynamics.  We  present  the  flow  equations  and  derive 
the  sensitivity  equations  as  for  the  1-D  model  problem.  An  adaptive  finite  element  tech¬ 
nique  is  described  and  used  to  solve  the  state  and  sensitivity  equations.  The  numerical 
approximations  are  presented  to  investigate  the  convergence  of  the  adaptive  grids. 


3.1  The  Navier-Stokes  Equations 


We  consider  steady-state  2-D  flow  of  an  incompressible,  viscous  fluid  in  a  bounded  domain 
D  C  with  a  Lipschitz-continuous  boundary  T.  Our  discussion  and  presentation  will 
follow  the  notation  in  [19]  with  the  exception  that,  in  most  cases,  vector  notation  will  be 
used  for  ease  of  presentation.  Let  u  =  [u(x,y),v(x,y)]T  represent  the  velocity  field,  and 
define  the  stress  tensor,  r,  by 

r(u)  =  y  [Vu  +  (Vu)7]  (3.1) 

where  y  denotes  the  viscosity.  Kundu,  in  [31],  gives  the  differential  form  of  the  principle  of 
conservation  of  mass  as 

dp 

+  V  •  (pn)  =  0,  (3.2) 

where  p  is  the  fluid  density.  A  fluid  is  incompressible  if  its  density  does  not  change  with 
pressure.  The  steady  incompressible  form  the  continuity  equation  then  becomes 


V  •  u  =  0. 


(3.3) 


24 
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The  steady-state  statement  of  conservation  of  momentum  is  the  Navier-Stokes  equation: 

-V  •  r(u)  +  VP  - 1-  p(u  -  Vu)  =  pf,  (3.4) 

where  P  denotes  the  ambient  pressure  and  f  represents  a  density  of  body  forces  per  unit 
mass  (e.g.  gravity).  In  general,  the  viscosity  p  is  a  function  of  the  temperature.  Since  we 
neglect  temperature  differences  for  the  problems  we  consider  herein,  p  is  assumed  to  be 
constant  and  can  be  taken  outside  the  derivative.  This,  along  with  the  fact  that  the  flow 
is  incompressible,  gives  the  following  identity: 

V  •  r(u)  =  pAu.  (3.5) 


As  commonly  done,  we  set 


(3.6) 


where  p  is  the  kinematic  pressure  and  v  is  the  kinematic  viscosity.  Then,  using  (3.5)  the 
system  (3.3)-(3.4)  becomes 


V  •  u  =  0  1 
— vAu.  +  Vp+  u  •  Vu  =f  J 


in  Cl. 


(3.7) 


It  is  useful  to  consider  a  non-dimensionalized  form  of  the  Navier-Stokes  equations.  To  do 
this,  define  a  length  scale  L,  a  velocity  scale  U,  and  a  reference  pressure  p0,  for  the  flow. 
The  dimensionless  variables  for  the  flow  can  be  defined  as  follows: 


X  =  Z’ 


u 


P  = 


P~Po 
pU2  ' 

The  non-dimensional  Navier-Stokes  equations  can  be  written  as 


u  =  — ,  and 


V  •  u 


pUL 


zlu  +  Vp  +  u  •  Vu 


in  Cl. 


(3.8) 


where  f  is  the  non  dimensionalized  f.  Dropping  the  tilde’s  and  defining  the  Reynolds 
number,  Re  =  we  have 


V  •  u  =  0  1  .  _ 

—jjfcAvi  +  Vp  +  u  •  Vu  =f  /  m  a  (3-9) 

For  ease  of  notation  and  to  maintain  consistency  with  [19],  replace  with  v  to  obtain  the 
working  version  of  the  Navier-Stokes  equation. 


V  •  u  =0  \ 
—vAu  +  Vp  -f  u  •  Vu  =  f .  J 


in  Cl. 


(3.10) 
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3.2  The  Homogeneous  Dirichlet  Problem 


We  begin  by  considering  the  case  of  the  homogeneous  Dirichlet  boundary  condition 

u  =  0  on  T.  (3.11) 

In  order  to  discuss  existence  and  uniqueness  and  to  introduce  the  variational  form  of  the 
problem,  we  introduce  some  standard  function  spaces  as  well  as  some  necessary  bilinear 
and  trilinear  forms. 


3.2.1  Function  Spaces  and  Notation 

As  usual,  we  let  L2(Q)  denote  the  space  of  square  integrable  functions  and  (-,  •)  and 
denote  the  L2  inner  product  and  norm,  respectively.  Let 


L2(n)  =  {qeL2(D):  JqdQ  =  0} 

Q 


The  Sobolev  spaces,  If1  and  Hq,  are 

ff'(Cl)  =  {v€  L*(n):  £  L2(fi)},  and 

Hq(CI)  =  {v  €  H1:  w|r  =  0}, 

with  inner  product,  (-,  •)],,  norm,  ||  -  ||l5  and  seminorm,  1*^,  defined  by 

/  \  /  \  ,  fiu  &v\  /9u  dv. 

(«,  n)i  =  (u,  v)  +  (—,  —)  +  — ) 


\u 


respectively. 


INK  =  (u,  u)\/2,  and 

I  —  ( t—  —\  i —  ® u\\l/ 2 
1  \  dx'  dx  dy’  dy  )  7 


Also,  denote  by  H  l(Q)  the  dual  space  consisting  of  bounded  linear  functionals  on  HUQ). 
The  norm  for  H~1(Q)  is  given  by 


IMI-i  =  sup 


(&  v) 


o^veH^n)  Mi 
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Recall  the  trace  theorem  (see  [10])  which  proves  the  existence  of  a  trace  operator  7  as 
follows: 

Theorem  1.  Let  Cl  be  bounded ,  and  suppose  Cl  has  a  piecewise  smooth  boundary.  In  addi¬ 
tion,  suppose  Cl  satisfies  the  cone  condition.  Then  there  exists  a  bounded  linear  mapping 

1 :  R1(SI) -H?( D,  IItMII  <  Mil,  (3.12) 

such  that  'jv  =  u|r  for  all  v  E  Cl(Cl). 


As  noted  in  [39],  the  range  of  the  trace  operator  is  H1^2(V).  The  norm  for  functions  g 
belonging  to  H1l2( f)  can  be  defined  by 


I  Ip  II 1/2  = 


inf 

veH^n), 
v—q  on  r 


IHIi. 


The  vector-valued  counterparts  of  these  spaces  in  R2  will  be  denoted  by  bold-face  symbols, 
i.e., 

H1^)  =  (Hl(Cl))2  =  {v  :  Vi  €  H\Cl)  for  *  =  1,2}  , 

H_1(0)  =  (H-'iCl))2  =  {v:Vi<=  H~l(Cl)  for  i  =  1,2}, 

H1/2(T)  =  (F1/2(r))2  =  {v  .  v.  e  Fi/2(r)  for  •  =  1?  2}  ,  etc. 

The  norm  for  Hx(fl)  is  defined  by 

/  2  \  1/2 

IMIi=(£NI? 

Vi=l 

The  divergence  free  subspace  of  H£(f2),  Z0,  is  given  by 

Z0  =  {veHx(O)|V-v  =  0}. 

We  define  the  following  bilinear  form: 

a0(u,v)  =  J i/Vu:Vvdf2  V  u,  veH1^) 
n 


£(§■£)■ 


where 


Dawn  L.  Stewart 


Chapter  3.  2-D  Flow  Problems 


28 


Also,  let 


and 


b(u,  q)  =  -  J  qV  •  udQ  Vue  H^ft),  V  q  e  L2(0), 

n 

ai(w;u,v)  =  J  (w  •  Vu)  -vdQ  V  u,  v,  w  e  Hx(ft). 


3.2.2  Existence  and  Uniqueness  of  Solutions  to  the  Variational 
Form 


We  present  some  basic  results  regarding  the  existence  and  uniqueness  of  solutions  to  the 
variational  form  of  the  homogeneous  Navier-Stokes  problem  as  presented  in  [19].  To  begin, 
we  note  that  the  homogeneous  partial  differential  equation  (3.10)-(3.11)  can  be  written  in 
variational  form  as  follows. 

Variational  Problem  1.  Given  f  e  H-X(f2),  find  a  pair  (u,p)  €  Z0  x  L%(n)  such  that: 

a0(u,v)  +  oi(u;u,v)  +  6(u,p)  =  (f,  v)  V  v  e  (3.13) 

Recall  that  the  trilinear  form,  gi(*;  -,  •)  has  some  nice  properties  described  in  the  following 
lemma. 

Lemma  1.  Let  u,  v  €  Hx(0)  and  let  w  e  Hx(fi)  with  V  •  w  =  0  and  w  •  n|r  =  0.  Then, 
the  trilinear  form  a1(-:  -,  •)  is  continuous  on  (H1( Q,))3  and  satisfies: 

ai(w;u,v)  +  oi(w;v,u)  =  0,  (3.14) 

ai(w;v,v)  =  0.  (3.15) 

Now,  the  existence  result  from  [19]  follows: 

Theorem  2.  For  f  e  H-1(fl),  there  exists  at  least  one  pair  (u ,p)  €  Z0  x  LZ(Q.)  which 
satisfies  (3.13). 

In  order  to  discuss  the  uniqueness  of  the  solutions  (u,p)  to  (3.13),  we  introduce  the  norm 
of  the  trilinear  form  oi(-;  •,  •),  denoted  N,  and  defined  by 

hr  ai(w;u,v) 

V  =  sup  t  .  .  .  ,  . 

u,v7w€Z0  \^\i  M-l  |w|x 
u,v,w^0 


(3.16) 
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We  also  set: 


llfll 


sup 

v€Zo 

v^O 


Theorem  3.  If  f  €  H-1(ft)  and 


Wh?)  llffl*.  <  l. 


then  the  Variational  Problem  1  has  a  unique  solution  fu,  p)  e  Z0  x  Z-o(fi). 


(3.17) 


(3.18) 


3.3  The  Nonhomogeneous  Dirichlet  Problem 


Now  consider  the  more  general  case  of  a  nonhomogeneous  Dirichlet  boundary  condition 

u  =  g  on  T.  (3.19) 

Denote  by  I\-,  0  <  i  <  p,  the  connected  components  of  the  boundary  T  as  depicted  in 
Figure  (3.1).  We  shall  henceforth  assume  that 


g  -  nds  =  0, 


0  <  i  <  p. 


(3.20) 
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The  variational  form  of  the  nonhomogeneous  partial  differential  equation  (3.10)-  (3.19)  is 
obtained  using  standard  weak  formulations.  In  particular,  we  consider  the  problem: 

Variational  Problem  2.  Given  f  €  H-1(f2),  find  a  pair  (u,p)  €  HX(Q)  x  L§(Q)  such 
that: 

G0  (u,  v)  +  «i  (u;  u,  v)  +  6(u,  p)  =  (f ,  v)  V  v  €  H£  (ft)  ’j 

V  •  u  =  0  in  Q,  i  (3.21) 

u  =  g  on  T.  J 

In  order  to  prove  existence  for  the  nonhomogeneous  problem,  we  need  the  following  tech¬ 
nical  result  due  to  Hopf  (see  Lemma  2.3,  page  287  in  [19]). 

Lemma  2.  Let  g  £  Hz  (F)  satisfy  (3.20).  For  any  e  >  0,  there  exists  a  function  u0  = 
uo(e)  £  Hx(f2)  such  that 

V  -  u0  =  0,  u0|r=g,  (3.22) 

|ai(v;uo,v)|  <  e|v|J  V  v  €  Z0  (3.23) 

The  following  existence  theorem  may  be  found  in  [19]. 

Theorem  4.  Let  f  £  H-1(f2)  and  g  €  Hz(r)  satisfying  (3.20).  There  exists  at  least  one 
pair  (u ,p)  £  Hx(0)  x  Z|(0)  which  is  a  solution  of  (3.21). 


Before  stating  the  uniqueness  result  again,  we  make  a  few  definitions.  For  any  function 
Uo  £  Hx(fi),  define 


and 


where 


p{  u0)  =  sup 


vez0 

v#0 


Qi(v;u0,v) 


IK(f;uo)||z'0  =  sup 
vez0 
v^O 


(3.24) 


(3.25) 


(/,  v)  =  (f,  v)  -  g0(u0,  v)  -  oi(uo;  u0,  v).  (3.26) 

Define  i/0  =  o0(Q,  f,  g)  as  in  [19]  by 

no  =  inf  |p(u0)  +  (^/’||i(f;  u0)||Z'o) " ;  u0  €  HX(S7)  satisfies  (3.22)| .  (3.27) 


The  basic  uniqueness  result  for  (3.21)  is  found  in  [19].  We  state  it  below  for  convenience. 
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Theorem  5.  Assume  the  hypothesis  of  Theorem  4.  If  u  >  v0(Q,  f,  g),  then  the  Variational 
Problem  2  has  a  unique  solution  (u ,p)  G  Hx(0)  x  Lq{Q). 

REMARK:  Although  the  results  above  provide  existence  and  uniqueness  for  the  basic 
Navier-Stokes  problems,  they  do  not  address  the  continuity  and  differentiability  of  these 
solutions  with  respect  to  the  parameter  q.  For  example,  in  the  problems  considered  below 
we  have  D  =  ft(q),  f  =  f(q),  and/or  g  =  g(q)  so  that  u0  =  i/0(fi,f,g)  = 
where  q  €  JRn  is  some  parameter  defining  the  flow.  We  need  to  establish  the  smoothness 
of  these  mappings  in  order  to  address  the  existence  and  uniqueness  of  solutions  to  the 
sensitivity  equations.  This  is  the  subject  of  the  following  sections. 


3.4  An  Abstract  Framework  for  Navier-Stokes 


In  this  section,  we  present  an  abstract  framework  for  analyzing  the  dependence  on  q  of 
solutions  to  the  nonhomogeneous  Dirichlet  problem  for  the  Navier-Stokes  equations.  We 
will  show  the  continuity  of  solutions  with  respect  to  parameters  for  two  specific  cases  and 
we  will  conclude  with  results  about  the  differentiability  of  those  solutions.  We  extend  the 
framework  in  [19]  to  certain  parameter  dependent  flows. 


3.4.1  The  Framework 

Let  X  and  X  be  two  Banach  spaces  and  Q  C  A  C  IRn ,  where  Q  is  open  and  A  is  compact. 
Given  a  Cp-mapping  (p  >  1) 

F  :  (q,u)  €  A  x  X  F(q,  u )  G  X,  (3.28) 

we  are  interested  in  solutions  to  the  state  equation 

F(q,u)  =  0.  (3.29) 

Let  {(g,  u(q));  q  G  A}  be  a  branch  of  solutions  of  equation  (3.29).  This  means  that 

q  ->■  u{q)  is  a  continuous  function  from  A  into  X ,  and  (3.30) 


F{q,u(q))  =  0. 


(3.31) 
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Moreover,  we  suppose  that  these  solutions  are  nonsingular  in  the  sense  that: 

DuF(q,  u(q))  is  an  isomorphism  from  X  onto  X,  Vq  £  A.  (3.32) 

As  an  immediate  consequence  of  (3.32),  it  follows  from  the  Implicit  Function  Theorem  (see 
e.g.  [40])  that  q  -4  u(q)  is  a  Cp-function  from  A  into  X. 

3.4.2  Using  the  Framework 

We  now  show  that  the  parameter-dependent  Dirichlet  problem  for  the  Navier-Stokes  equa¬ 
tions  in  the  velocity-pressure  formulation  (3.10)  -  (3.19)  fits  into  this  abstract  framework. 
We  assume  that  any /all  of  the  following  hold: 

g  =  g(q),  and/or  (3.33) 

f  =  f(q)  (3.34) 

where  q€QC4C|fisa  design  parameter  for  the  flow.  Define 

X  =  X  =  Hx(fl)  x  L\{D),  (3.35) 

and  the  intermediate  space 


Y  =  H~1(Q)xV,  (3.36) 

where  V  —  |g  €  H1/2(r);  g  •  n  ds  =  0,0  <  i  <  p|.  Next  we  define  a  linear  operator  T 

as  follows:  given  (f„g*)  6  7,  we  denote  by  (u *,/>*)  =  T(f*,g,)  €  X  the  solution  of  the 
Dirichlet  problem  for  the  Stokes  equations: 


V  -  u*  =  0  in  D  \ 

— Zlu*  -I-  Vp*  =  f,  in  Q  >  . 

=  g*  on  T  J 

In  addition,  let  P  :  Q  -4  Y  be  defined  by 

p(q)  =  (^f(q)>g(q)) , 

and  the  non-linear  operator  MC  :  X  -4-  Y  be  given  by 

= (v,y))  =  (l  (v,fi + V2ij)  >  °)  ’ 


(3.37) 


(3.38) 


(3.39) 
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where  v  is  the  constant  ~  as  before. 

Now  finally,  with  the  data  (f ,  g)  €  Y  we  associate  a  ^-mapping  G  from  Qx  X  into  Y 
defined  by 


G  :  (q,  v)  -4  G( q,  v)  =  NC(v)  -  P( q)  (3.40) 

and  we  set 


F(q,v)  =  v +  TG(q,v).  (3.41) 

The  following  result  follows  directly  from  Lemma  3.1  in  [19]. 

Lemma  3.  The  pair  (u(q),p(q))  €  Hx(f2)  x  Lg(fi)  is  a  solution  of  (3.10)  -  (3.19)  if  and 
only  if  (q,u(q))  is  a  solution  of  (3.29)  where  u(q)  =  (u(q),p(q)/i/)  and  where  the  spaces 
X  and  X  are  defined  by  (3.35)  and  the  compound  mapping,  F,  is  defined  by  (3.41). 


3.4.3  Continuity  of  Solutions  with  Respect  to  Data 

We  now  address  the  continuity  of  solutions  (u(q),p(q))  to  (3.10)-(3.19)  with  respect  to 
changes  in  the  parameter  q.  We  assume  the  map  q  e  Q  ->•  (g(q),f(q))  €  V  x  H-1(Q) 

is  continuously  Frechet  differentiable.  Note,  that  this  is  certainly  true  for  the  cylinder 

problem  presented  in  §3.7.1. 

To  begin,  we  need  Lemma  1.3.2  from  [19]. 

Lemma  4.  There  exists  a  continuous  linear  function  V  :  V  -»  H^fl)  such  that  for  each 
g  €  V,  we  have  w  =  X>(g)  satisfies 

V  -  w  =  0,  and  wl  =  g  )  , 

IMIi  <  C\\z\\i/2  }‘  (3‘42^ 

The  following  corollary  is  a  direct  consequence  of  Lemma  4. 

Corollary  1.  The  map  from  q  €  Q  2>(g(q))  €  Hx(ft)  is  Frechet  differentiable. 


In  order  to  analyze  the  parameter-dependent  solution  to  the  weak  form  of  the  nonhomo- 
geneous  Navier-Stokes  problem,  we  need  a  result  analogous  to  Lemma  2  for  parameter- 
dependent  boundary  functions.  The  following  result  may  be  established  by  a  straight 
forward  extension  of  Lemma  2.3  in  [19]. 
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Lemma  5.  There  exists  a  continuous  linear  function  T :  A  x  ]R+  ->  Hx(ft)  such  that  for 
each  q  G  A,  e  >  0 

^(q,e)  =  u0(q)  (3.43) 

where  u0(q)  is  defined  by  Lemma  2  and  satisfies 

V  •  Uo(q)  =  0,  u0(q)|r  =  g(q),  (3.44) 

and 

Mv;uo(<l)»v)|<e|v|J  (3.45) 

for  each  v  G  Z0. 

We  now  consider  continuity  with  respect  to  the  right  hand  side  function  f.  We  again 
consider  an_abstract  framework  for  the  Navier-Stokes  equations.  We  define  a  map  Q  : 
X  xy  ->  Z.  Here,  X  is  the  set  of  forcing  functions,  f  G  H_1(ft),  and  y  is  defined  as 
follows: 

y  =  { (u, p) :  u  €  Hj(ft),  and  p  G  H1{D)/1R}  .  (3.46) 

The  range  of  Q  is  contained  in 

Z  =  H-1(fi)  x  Ll(Q)  (3.47) 

and  Q  is  defined  by 

Q(x  =  f ,  y  =  (u, p))  =  {-uAn  +  (u  •  V)u  +  Vp  -  f ,  V  •  u) .  (3.48) 

Note  that  Q  G  C1  and  the  Frechet  derivative  at  a  point  (x0  =  f0,  Vo  =  (u0,p0))  is  given  by 

IDQ(X 0, «,)]  (*,  y)  =  (  ~f +  <Uo  •  V>u  +  ;J)uo  +  Vp  -  |  (3  49) 

It  is  also  clear  that  (u,p)  satisfies  the  homogeneous  Navier-Stokes  equations,  (3.10)  -  (3.11), 
with  right  hand  side  f  if  and  only  if  Q  (x  =  f,  y  =  (u,  p))  =  0. 

The  Implicit  Function  Theorem  implies  that  (u,p)  is  a  continuously  differentiable  function 
off  if  the  linear  map  \DyG{xQ,  y0)}  is  an  isomorphism.  But  this  is  equivalent  to  the  condition 
that  the  homogeneous  sensitivity  equation  has  a  unique  solution  in  y  for  each  f  G  if-1  (ft). 
Weshow  in  the  following  that  the  sensitivity  equation  does  indeed  have  a  unique  solution 
in  3^,  so  that  we  have  (u,p)  is  a  C 1  function  of  f.  We  will  then  return  to  the  smoothness  of 
solutions  to  the  parameter-dependent  Navier-Stokes  equations. 
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3.5  Analysis  of  the  Sensitivity  Equations 

We  begin  by  stating  a  general  form  of  the  sensitivity  equations  for  the  parameter-dependent 
Navier-Stokes  equations.  We  show  that  if  we  have  a  unique  solution  (u ,p)  of  the  Navier- 
Stokes  problem,  then  we  have  a  unique  solution  (s,  r)  of  our  sensitivity  equation.  Lastly, 
we  will  use  an  abstract  formulation  of  the  Navier-Stokes  problem  and  the  implicit  function 
theorem  to  show  that  the  solution  (u ,p)  is  in  fact  a  nonsingular  solution  of  the  Navier- 
Stokes  problem. 


3.5.1  A  General  Formulation  of  the  Sensitivity  Equations 

Let  q  €  2R+  represent  some  fixed  sensitivity  parameter  q{,  (1  <  i  <  N).  Now  denote  the 
sensitivity  forcing  function,  J|( q ),  by  f*  and  the  boundary  function,  J J($),  by  g*.  Lastly, 
recall  that  in  some  cases  Q  =  Cl(q).  Then  as  in  the  case  of  the  1-D  model  problem,  the 
sensitivity  equations  are  obtained  by  implicitely  differentiating  the  flow  equations  and  their 
associated  boundary  conditions.  The  sensitivity  equations  fit  into  the  following  general 
form.  Given  f*  €  H-1(ft)  and  g*  6  H>(T)  satisfying  (3.20)  and  (u ,p)  a  solution  of  (3.10)- 
(3.19),  find  a  pair  (s,  r)  €  HX(Q)  x  Lg(f2)  such  that  (s,  r)  satisfies 

V-s  =0  in  Cl  1 

-vAs  +  Vr  +  u  -  Vs  +  s  -  Vu  =  f *  in  Cl  >  .  (3.50) 

s  =g*  on  T  j 

We  can  write  the  variational  form  of  (3.50)  as: 

Variational  Problem  3.  Given  f*q  €  H_1(0)  and  g*  €  H^(T)  satisfying  (3.20)  and  (u ,p) 
a  solution  (3. 10)-(3. 19),  find  a  pair  (s,r)  6H'(fi)  x  L2q(Q)  such  that: 

a0(s,z)  +ai(u;s,z)  +a1(s;u,z)  +  6(z,r)  =(f*,z)  'j 

V  •  s  =  0  in  Cl  >  ,  (3.51) 

s  =g*  on  T  j 

Vz  €  Hj(fl). 


3.5.2  Existence  and  Uniqueness  of  Solutions  to  the  Sensitivity 
Equations 

We  state  and  prove  the  following  existence  and  uniqueness  result  following  [19]. 
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Theorem  6.  Assume  the  hypotheses  of  Theorem  5.  If  (u..  p)  is  a  unique  solution  of  the 
Variational  Problem  2,  then  there  exists  a  unique  solution,  (s  ,r),  to  Variational  Problem 
3. 

Proof.  In  order  to  check  that  the  equations  (3.51)  have  a  unique  solution,  it  is  sufficient  to 
prove  that  the  bilinear  form 

(s,  z)  -4  c(s,  z)  =  a0 (s,  z)  +  Oi  (u;  s,  z)  +  ax (s;  u,  z)  (3.52) 

is  V-elliptic.  But  it  follows  from  (3.15)  that 

c(s,s)  =  i'|s|i  +  «i(s;u,s)  VsgZo-  (3.53) 

Now,  assume  that  v  >  i/Q  where  uq  is  defined  by  (3.27).  Then,  there  exists  a  function 
u0  €  H1(fl)  such  that 


V-uo  =  0,  u0|r  =  g* 
v>p(uo)  +  (^||Z(r;u0)||z{()1. 
Next,  setting  u  =  u0  +  w,  we  have  by  (3.16)  and  (3.24) 

K(s;u,  s)|  <  |ai(s;u0,s)|  +  |oi(s;  w,s)| 
<  (p(u0)  +A/’|w|1)  |s|j  . 

Since 

[K(f;no)lk 


we  obtain 


Hi  < 


c(s 


,  s)  > 


p(u o)  - 


JW;«o)ll* 


v  -  p(u0) 


N? 


=  w 

so  that  the  ellipticity  property  holds. 


^•n1-  k-,(uo))2  )|s|1 


(3.54) 


3.6  Differentiability  of  Solutions  with  Respect  to  g 


We  return  now  to  the  smoothness  of  solutions  to  the  parameter-dependent  Navier-Stokes 
equations. 
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Note  that  the  parameter  i/0  =  z/0(q)  is  a  continuous  function  of  q  in  this  case  since  the 
map  from  q  -4  u0  is  continuous  and  the  map  from  u0  -4  p  is  continuous.  With  this,  the 
continuity  of  the  map  (f,  g)  -4  (u,p),  and  the  fact  that  q  -4  (f(q),  g(q))  is  C\  we  can  now 
establish  the  following  result. 

Theorem  7.  Assume  that  q  £  Q  is  such  that^u  >  z'o(q)-  There  exists  a  neighborhood ,  Q 

such  that  for  all  q  €  Q  the  solution  of  the  variational ,  parameter- 
dependent  Navier-Stokes  equations  (3.21),  with  f  =  f  (q)  and  g  =  g(q)  exists  and  is  unique. 
Moreover,  the  solution  is  a  C1  function  of  q. 


Proof:  We  consider  the  map  F  defined  by  (3.41).  We  have  shown  that  there  exists  Q  C  TP+ 
so  that  we  have  a  branch  of  solutions  {(q,u(q));q  €  q},  i.e.  that  both  (3.30)  and  (3.31) 
hold  for  q  €  Q,  for  both  cases  presented  in  §3.4.3  above. 


We  now  turn  our  attention  to  the  Frechet  derivative  of  F,  DUF.  We  have 

[DuF(q,u)](q,u)  =  u  +  T  [DuG(q, «)]  (q,  u) 

=  u  +  T[DnMl{q,u)}{u) 


=  u  +  T 


(3.55) 

(3.56) 

(3.57) 


Again,  the  fact  that  DUF  is  an  isomorphism  can  be  shown  to  be  equivalent  to  the  fact  that 
the  homogeneous,  sensitivity  equation  has  a  unique  solution  for  all  q  G  Q.  We  have  shown 
in  §3.5.2  that  the  sensitivity  equation  does  indeed  have  a  unique  solution.  Then,  by  the 
Implicit  Function  Theorem,  we  have  that  q  -4  (u(q),p(q))  is  a  C1  function  from  Q  into  X. 


3.7  Two  Specific  Problems  and  Their  Sensitivity 
Equations 


We  now  return  to  analyzing  design  sensitivities  and  focus  on  two  specific  parameter- 
dependent  flow  problems  along  with  their  sensitivity  equations. 


3.7.1  Flow  around  a  Cylinder 

We  begin  by  considering  the  standard  problem  of  2-D  flow  around  a  cylinder.  This  prob¬ 
lem,  a  non-homogeneous  Dirichlet  problem  on  a  bounded  domain,  is  modeled  on  the  prob- 
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(L,l) 


(U-l) 


lem  where  the  boundary  is  an  infinite  strip.  We  assume  parabolic  inflow  into  a  chan¬ 
nel  containing  a  cylindrical  obstruction  whose  geometry  is  shown  in  Figure  (3.2).  Here, 
^  =  [— 2,  L\  x  [—1,1],  where  L  €  1R+  is  a  fixed  number.  We  assume  that  L  is  large  enough 
for  the  outflow  to  have  returned  to  the  same  parabolic  velocity  profile  present  at  the  inflow. 
The  governing  equations  are  the  2-D  incompressible  Navier-Stokes  equations  presented  ear¬ 
lier  (see  (3.10))  with  f  =  0.  The  boundary  conditions  at  the  inflow  and  outflow  are  given 
by 


u(~2,  y ;  q )  =  u (L,  y,  q )  = 


(i  +  ?)(i  -y2) ' 
0 


(3.58) 


for  —1  <  y  <  1,  where  q  £  R+  is  a  parameter  describing  the  strength  of  the  inflow.  No¬ 
penetration  and  no-slip  conditions  are  applied  on  the  top,  bottom,  and  cylinder  sidewalls 
(i.e.  u  =  0). 


We  are  interested  in  calculating  sensitivities  with  respect  to  the  inflow  parameter  q.  As 
before  with  the  case  of  the  1-D  model  problem,  define  the  sensitivity  as  follows: 

— w(a:,  y,  q),  —v(x,  y,  q)  \  =  s(x,  y ;  q)  =  s(z,  y).  (3.59) 

The  sensitivity  equations  are  obtained  by  implicitly  differentiating  the  system  (3.10)  and 
its  associated  boundary  conditions  with  respect  to  the  parameter  q.  Assuming  the  order 
of  differentiation  can  be  interchanged,  we  obtain  the  following: 

Vs  =  0  | 

+  V  +  u  •  Vs  +  s  •  Vu  =  0  J  in  a 
Observe  that  this  sensitivity  equation  does  not  assume  a  fully  developed  flow. 


(3.60) 
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m 


Figure  3.3:  Geometry  for  Flow  over  a  Bump 

The  sensitivity  boundary  conditions  at  the  inflow  and  outflow  are  obtained  by  differentiat¬ 
ing  (3.58)  producing 

*(-2,y;q)  =  s(L,y;q)=  ^  “ V  ^  ,  -1  <  y  <  1.  (3.61) 

After  differentiation,  the  no-penetration,  no-slip  conditions  for  u  imply  no-penetration,  no¬ 
slip  conditions  for  s.  In  §3.5  ,  we  present  the  weak  form  of  these  sensitivity  equations  and 
prove  an  existence  result. 


3.7.2  Flow  over  a  Bump 

We  also  consider  a  problem  examined  by  Burkardt  in  [11],  In  particular,  the  problem  is  2-D 
incompressible  flow  over  a  bump  in  a  channel.  The  geometry  of  the  channel  is  indicated 
in  Figure  (3.3)  with  Cl  =  [0,  L]  x  [0,3],  where  L  ^  0  is  again  a  fixed  number  representing 
the  length  of  the  channel.  As  before,  the  governing  equations  for  this  problem  are  the 
2-D  incompressible  Navier-Stokes  equations  presented  earlier  (see  (3.10))  with  f  =  0.  The 
boundary  conditions  for  the  flow  are  as  follows: 

u(0,  y)  =  u (L,  y)  =  J  ~V)  J 


(3.62) 
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for  0  <  y  <  3,  where  A  =  0.5  was  a  constant  parameter  describing  the  strength  of  the 
inflow.  Again,  no-penetration  and  no-slip  conditions  are  applied  on  the  top  and  bottom 
channel  sidewalls  as  well  as  on  the  bump. 


For  this  application,  we  examine  a  shape  sensitivity.  Here,  the  shape  of  the  bump  is  a 
cubic  spline  determined  by  a  parameter,  q.  We  examine  two  cases:  q  =  qx  €  R1  and 
q  =  (hi  hi  h)  €  JFP ■  We  seek  to  find  the  sensitivity  of  the  flow  in  the  channel  to  changes 
in  q,  which  we  denote  *  =  =  («„,»„). 

We  have  to  solve  a  separate  set  of  linear  sensitivity  equations  for  each  sensitivity.  The 
sensitivity  equations  are: 

V-Si  =0  ] 

-uAsi  +  V  (j^pj  +  u  •  Vsi  +  Si  -  Vu  =  0  J  m  (3‘63) 


In  this  case,  the  sensitivity  boundary  conditions  are  given  by 


uqi(x,  y) 

»«(*>  v) 


uqi(x,  y) 


vgi(x,  y) 


0  at  the  inflow,  outflow,  top,  and  bottom  walls; 
0  at  the  inflow,  outflow,  top,  and  bottom  walls; 
du  dh(x,  q) 

dy  % 

dv  dh(x ,  q) 


dy  % 


on  the  bump; 
on  the  bump; 


(3.64) 

(3.65) 

(3.66) 

(3.67) 


where  h(x,  q)  denotes  the  ^-coordinate  of  the  boundary  of  Q  for  1  <  x  <  3.  We  generate 
h{x,q)  using  a  cubic  spline  with  free  boundary  conditions.  If  q  e  R1,  qx  specifies  the 
height  of  the  spline  at  x  =  2.0.  In  this  case,  h(x,  q)  satisfies  the  following: 


h(x ,  q)  is  piecewise  cubic, 

Mi,  q)  =  M3,  q)  =  o, 

h(  2,  q)  =  qi 

h'(l,q)  =  h'(2,q)  =  h'(3,q)=0 


(3.68) 


For  q  €  JR3,  qu  q2,  and  qz  specify  the  height  of  the  spline  at  x  =  1.5, 2.0,  and  2.5,  respec¬ 
tively.  Here,  h(x,  q)  satisfies 


h(x,  q)  is  piecewise  cubic,  ' 

h(l,q)  =  h{3,  q)  =  0, 
h(1.5,q)  =  qx 

h( 2-0,  q)  =  q2  > . 

h(  2.5,  q)  =  qz 
h'  (1,  q)  =  h'  (3,  q)  =  0 

M(x,q)  is  continuous  at  x  =  1.5, 2.0,  and  2.5 


(3.69) 
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3.8  A  Finite  Element  Formulation 


Since  exact  solutions  to  the  Navier-Stokes  equations  are  unavailable  except  for  the  simplest 
of  problems,  we  now  turn  our  attention  to  finding  good  methods  for  computing  approxi¬ 
mations  to  the  solutions  of  the  flow  equations.  We  will  use  a  finite  element  approach  to 
solve  our  variational  problems  and,  later,  our  weak  sensitivity  equations  as  well. 


3.8.1  The  Discrete  Variational  Problem  and  Finite  Element 
Spaces 


To  begin,  we  choose  a  family  of  finite  element  spaces,  VA(Q)  and  SA(Q),  for  approximating 
velocity  and  pressure,  respectively.  Here  h  is  some  measure  of  the  size  of  the  grid  used  to 
divide  £2  into  finite  elements.  Then  the  finite  dimensional  weak  form  of  the  homogeneous 
Dirichlet  problem  becomes: 

Variational  Problem  4.  Given  f  €  H_1(£2),  find  a  pair  (u A,pA)  e  Vj}(S2)  x  Sft(Q)  such 
that: 

ao( uA,  vA)  +  fli(uA;  uA,  vA)  +  b(vA,pA)  =  (f,  vA )  V  vh  £  Vft(Q)  1 

b( uA,gA)  =  0  V  qA  £  SjfiD)  f  (3-70) 


When  solving  Navier-Stokes,  the  inclusions  of  our  finite  element  spaces,  V£  C  Hj(fi)  and 
S&  C  lg(0)  in  our  infinite  dimensional  function  spaces  are  not  by  themselves  sufficient  to 
produce  a  stable,  meaningful  approximation.  For  stability,  the  finite  element  spaces  must 
also  satisfy  the  additional  div-stability  or  inf-sup  condition.  Three  equivalent  forms  of  this 
condition  from  [22]  are  given  below. 

A  finite  element  space  is  said  to  satisfy  the  div-stability  or  inf-sup  condition  if  and  only  if 
one  of  the  following  equivalent  conditions  holds: 


•  Given  any  qA  £  S(fi 


sup 

o#v*ev§ 


V  |v*|,  )- 


(3.71) 


where  the  constant  7  >  0  may  be  chosen  independent  of  h  and  of  the  particular  choice 
of  j*  e 
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•  Given  any  qh  £  Sft  there  exists  vA  €  Vq  such  that 

KVW)  >7llff*llo  Iv^Ij,  (3.72) 

where  the  constant  7  >  0  may  be  chosen  independent  of  h  and  of  the  particular  choice 
of  qh  £  S%. 


There  exists  a  7  >  0,  independent  of  h,  such  that 

(  KvhAh) 


inf  sup  ■„  v.  ,  ’V  >  7. 
*&&&  <&**&/&  V Iloilo  |VA|  J 


(3.73) 


As  noted  in  [22],  the  div-stability  condition  (3.71)  -  (3.73)  loosely  ensures  that  as  h  ->•  0 
discretely  divergence  free  functions  tend  to  divergence  free  functions. 

The  discrete  Weak  Formulation  4  is  solved  in  the  primitive  variables,  uh  and  ph.  The 
weak  form  is  discretized  using  the  7  node  Crouzier-Raviart  triangular  element  (see  Figure 
(3.4))  which  is  a  type  of  “bubble”  element  described  in  [14].  This  element  uses  a  so-called 
enriched  quadratic  velocity  interpolant  and  a  discontinuous  linear  pressure.  Then,  if  %  is 
some  triangulation  of  our  domain,  ft,  the  finite  element  space,  Sq  can  be  formally  defined 
as  follows: 

So  =  I?  :  Q  €  Si  (A),  A  €  7ft;  J  q  dft  =  0  j  ,  where  (3.74) 

where  Sx( A)  is  the  space  of  linear  polynomials  on  the  triangle,  A.  The  velocity  finite 
element  space  is  the  space  of  continuous  piecewise  quadratic  polynomials  augmented  in 
each  triangle,  A  €  7ft,  by  the  cubic  bubble  function  that  vanishes  on  the  three  edges  of  A. 
Thus, 

Vo  =  {v:v€[S2(A)©B2(A)],A€7I;  v  €  C0(ft);  v  =  0onr},  (3.75) 
where,  for  any  A  £  7ft, 

B2(A)  =  {v  £  82(A)  :  v  =  A1A2A3U,  u  £  So(A)}  .  (3.76) 

Here,  A*(x),  i  =  1,2, 3,  denote  the  barycentric  coordinates  of  the  point  x  £  A  with  respect 
to  the  vertices  of  A,  see  [14]. 


3.8.2  The  Penalty  Method 

The  pressure  degrees  of  freedom  are  treated  by  the  penalty  method  as  in  [17].  The  penalty 
method  allows  for  the  uncoupling  of  the  momentum  and  continuity  equations  through  the 
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X  Pressure  Node 

(1  nodal  point  +  2  derivatives) 

O  Velocity  Node 
(7  nodal  points) 


Figure  3.4:  Crouzier-Raviart  element 


use  of  a  small  perturbation  parameter,  e.  This  provides  large  computational  advantages 
as  it  reduces  the  size  of  the  system  of  equations  being  solved  and  eliminates  the  need  for 
partial  pivoting  when  solving  the  resultant  linear  systems. 

The  idea  of  the  penalty  method  is  to  perturb  the  continuity  equation  with  a  small  term 
containing  the  pressure: 


sp  +  V  -  u  =  0.  (3.77) 

One  can  consider  this  perturbation  as  the  introduction  of  a  slight  artificial  compressibility. 
With  this  perturbation,  the  pressure  p  is  expressed  as 


p  =  -/?V  •  u,  (3.78) 

so  that  once  a  solution  for  u  is  known,  p  is  determined  by  the  relation  (3.78).  Here  /3  =  1/e 
is  called  the  penalty  parameter.  With  (3.78),  we  can  eliminate  pressure  from  the  discrete 
form  of  the  momentum  equation.  It  can  be  shown  that  a  solution  (u/3,^)  of  the  system, 

— /?V  -u  =p  1  _ 

—vAvl  +  Vp  +  u  •  Vu  =  f  J  m  (3-79) 

approaches  the  solution  of  the  unperturbed  system  (3.10)  as  p  ->  oo  (see  [37]  for  details 
for  the  homogeneous  and  nonhomogeneous  cases). 
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3.8.3  An  Implementation  Method 


Let  <^(x),  j  =  1,2,..., TV  and  ^(x),  l  =  1,2,...M  be  basis  functions  for  the  finite 
element  spaces  Vq  and  respectively.  The  Galerkin  approximations  (u h,ph)  of  the 
penalized  variational  form  can  be  written  as 


uA(x)  =  ]CUA(X)>  or,  alternately  f 

j=i  \v  W  J  \ 

M 

ph(*)  = 


(3.80) 

(3.81) 


1=1 


Substituting  these  approximations  into  the  penalized,  discrete  weak  form,  we  obtain  the 
following  system  of  Galerkin  equations: 

aQ(uh,  fy)  +  ax(uA;  uft,  fy)  +  ph)  =  (f,  <f>.)  j  =  1, 2, . . . ,  N  and  (3.82) 


/%(u\</>/)  =  (ph,  ipi)  /  =  1,2, 


(3.83) 


Using  the  definitions  of  ao,  ax,  and  b  and  separating  the  momentum  equations  into  u,  v 
components  for  clarity,  we  can  write  (3.82)  -  (3.83)  in  integral  form  as: 


J  vVuh  ■  V<f>i  + 


duh 

h,duh\ 

dQ  = 

f 

dx 

~ty)' 

4>i-ph^r 

ox 

/  /ifc. 
Jn 

,dvh 

hdvh\ 

J  hd<t>i 1 

dO  = 

f 

dx 

~dy)' 

<t>i-Ph~zr 
dy  J 

/  /2&< 
Jn 

h<t>i<m,  *  =  1,2,..., IV  (3.84) 


'  duh 

dvh' 

dx 

+  ~dy. 

iptdCl 


=  f  p%<m ,  /  =  1, 2, . . . ,  M.  (3.86) 

J 17 


Equation  (3.86)  can  be  written  more  succinctly  in  matrix-vector  notation  as 

-f3LTuh  =  Dph  (3.87) 

where, 

ph  €  Rm ,  ph  =  [pi,p2, . . .  ,pMf ,  (3.88) 

nheJR2N,  uh  =  [ui,u2,...,un,vi,v2,...,vn]t  ,  (3.89) 
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and  D  e  ]RMxM  is  the  pressure  mass  matrix 

Dij  =  /  dD  i,  j  =  1, 2, . . . ,  M, 

Jo 

and  L  €  1 RMx2N  is  the  continuity  matrix  defined  as  follows: 

L..  =  I  i  =  1,2, . . . ,  M,  j  =  l,2,...,N 

V  {fn&dTdn  i  =  l,2,...,M,  j  =  N+l,N+2,. 

We  can  solve  (3.87)  for  so  that 

ph  =  -f3D~lLT\ih. 


(3.90) 


. . ,  2N 


(3.91) 


(3.92) 


To  linearize  (3.84)  -  (3.85),  we  consider  an  iterative  method  where  we  write  the  solution  at 
the  new  level,  n,  as  the  sum  of  the  preceding  level,  n  -  1,  and  a  correction: 


Uh,n  =  uA’ra_1  + 


(3.93) 


Substituting  (3.93)  into  the  convective  terms  of  (3.84)  -  (3.85)  and  neglecting  the  quadratic 
terms  in  du,  the  linearized  form  of  the  convective  terms  of  the  discrete  momentum  equations 
become 

r  + 

J n  L  dx  dx  dx 

v  sr + v  ~»r  ~  ”  “sH  ****  (3-94) 

/  + 

Jo  [  dx  dx  dx 

v  ~dT  +  v  ~W  ~ v  (3’95) 


for  i  =  1, 2, . . . ,  N. 

In  matrix-vector  form,  the  linearized  convective  terms  can  be  written  as 

N(uh,n~1)uh,n  -  g, 

where  N(u h’n~1)  €  JR2Nx2N.  The  matrix  N  is  defined  as 


(3.96) 


Nuu 

:  N 

■  1  yuv 

Nvu 

:  N 

•  1  yvv 

(3.97) 
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with  the  submatrices  defined  by  the  following: 

=  l  + “‘-'f + )  *  «• 

dvk,n 


lw"lw  ~L{ 


—  f 

J  Jit 


<f>: 


duh,n~  i 


Qvh,n-1 


g  = 


91 

92 


with 


(3.98) 

(3.99) 

(3.100) 

N 

9z  =  L  («*,n_lfl!15ri  +  wM"lfi2Szi  j  <t>i  dn  i  =  1, 2, . . . ,  N. 

(3.101) 


4>i  d£l,  and  [IV**]^-  =  f  $/- 


'<Pi  dQ  5 


J  dy  l-TOJv  JM3  dx 

for  i,j  =  1, 2, . . . ,  N.  The  vector  g  €  ]R2N  is  given  by 

9i  =  /„  („M-i4>£i  +  *  <iO  i  =  1, 2 . 


With  these  definitions,  the  complete,  linearized  momentum  equations  (3.84)  -  (3.85)  can 
be  written  in  matrix-vector  notation  as 


Suh’n  +  lV(uft’n_1)u'1’”  +  tLD~1Lt  uft>n  =  h  +  g, 
where  S  6  lR2*Kx~N  is  the  diffusion  matrix  with 

r  i//0v*.v^  =  1,2, ... ,iv 

Sy  =  \  v In  V(f>i_N  .  V^-_jv  i,j  =  1V+1, 1V+2, . .  .,21V 
(.  0  otherwise. 

and  the  right-hand  side  vector  h  €  M2N  defined  by 

hi 


(3.102) 


(3.103) 


h  = 


ho 


with  /  ^  /f 

l  =  Jn 


=  Jn Mi  dn  ?  =  i, 2, ... ,iv 

MidQ  i  =  1,2,...,  IV. 


(3.104) 


We  finish  this  section  with  a  few  notes  about  the  structure  of  the  matrices  L,  5,  and  D  for 
the  Crouzier-Raviart  element  described  in  §  3.8.1  above.  First,  note  that  with  the  use  of 
the  discretized  continuity  and  momentum  equations  and  Gauss’  theorem  that  the  velocity, 
and  in  turn,  the  pressure  derivatives  at  the  centroid  of  the  element  can  be  written  in  terms 
of  the  other  nodal  quantities.  This  results  in  a  velocity  approximation,  uft,  determined  by 
six  nodal  points  and  a  pressure  approximation,  ph,  determined  by  one  nodal  point.  The 
numerical  advantage  of  this  element  then  is  that  the  number  of  equations  and  unknowns 
in  each  element  is  reduced  by  3,  without  affecting  the  accuracy  of  the  approximation.  In 
addition,  the  matrix  D  can  be  written  as  a  diagonal  matrix  with  entries  D^  representing 
the  area  of  element  e*.  For  further  details  on  the  implementation  of  the  penalty  function 
method  with  the  Crouzier-Raviart  element,  the  reader  is  referred  to  [17]. 
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3.8.4  Adaptive  Methodology 


The  basic  idea  of  adaptive  gridding  is  to  use  some  error  estimation  technique  to  evaluate 
the  quality  of  your  finite  element  approximation  and  to  strategically  modify  the  grid  based 
on  that  evaluation.  The  grid  modification  scheme  allows  the  user  control  over  element  size 
and  grading.  This  process  has  been  shown  to  be  successful  in  resolving  shear,  stagnation 
points,  jets,  and  wakes  (see  [34],  [33],  [24],  [29],  [25]).  The  two  main  elements  of  the  adaptive 
process  are  error  estimation  and  grid  generation.  We  discuss  each  of  these  below. 


Error  Estimation 


The  error  estimation  is  performed  using  an  approach  introduced  by  Zhu  and  Zienkiewicz 
(see  [41],  [42],  [1])  and  involves  the  postprocessing  of  stresses  and  strains.  Recall  that  the 
energy  norm  of  u  is 


Ml* 


or,  given  in  Cartesian  coordinates, 


(3.105) 


Note  that  the  energy  norm  has  a  very  similar  form  to  the  H 1  seminorm.  In  fact,  it  can  be 
shown  that  they  are  equivalent  norms.  Using  the  energy  norm  for  the  velocity,  define  the 
so-called  Stokes  norm  of  the  solution  as 

ll(u,p)||s  =  \/lM|!  +  ||p||2.  (3.107) 

As  pointed  out  in  [35]  the  use  of  the  energy  norm  over  the  Hl  seminorm  offers  some 
advantages  especially  to  the  engineering  community.  Note  that  both  the  velocity  and 
pressure  norms  are  expressed  in  terms  of  surface  forces  which  are  the  quantities  of  prime 
interest  in  engineering  fluid  mechanics.  Secondly,  errors  computed  in  these  norms  can  be 
interpreted  as  errors  in  the  stresses  which  can  then  easily  be  related  to  errors  in  global 
quantities  such  as  lift  and  drag. 

The  Zhu  and  Zienkiewicz  approach  uses  the  Stokes  norm  to  measure  the  error,  e(u ,p)  = 
u  ,  Pex  P  ) ;  where  (u exiPex)  is  the  exact  solution  of  the  flow  problem  and  (u^,  p^1) 
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is  the  finite  element  approximation.  Then  the  norm  of  the  error  is 

l|e(u,p)||5  =  ^/l|u«-uA|||  +  ||pes-pA||2. 


(3.108) 


We  concentrate  on  forming  an  approximation  to  ||u«*  -  uA||£.  Since  the  exact  solution, 
and  more  particularly  the  gradients  of  the  exact  solution  are  not  known,  the  approach  is  to 
use  the  finite  element  approximation  to  construct  approximations  to  these  gradients.  Note 
that  the  finite  element  approximations  to  the  gradients  are  discontinuous  across  element 
faces  while  the  exact  gradients  are,  in  most  cases,  continuous  across  the  domain.  Thus, 
the  first  goal  in  error  estimation  is  to  obtain  continuous  approximations  to  the  discontin¬ 
uous  finite  element  gradients.  Two  methods  have  been  evaluated  for  this  process:  global 
projections  and  local  least  squares  projections.  Global  projections  are  done  over  the  entire 
domain,  fl,  and  involve  finding  the  best  approximation  to  the  discontinuous  finite  element 
gradients  in  the  original  continuous  finite  element  space.  For  example,  if  a  piecewise  linear 
approximation  of  the  flow  solution  is  calculated,  then  the  finite  element  gradient  is  piece- 
wise  constant  and  discontinuous  across  elements.  The  global  projection  would  replace  the 
piecewise  constant  gradient  approximation  with  a  piecewise  linear  approximation  calcu¬ 
lated  by  projecting  the  piecewise  constant  function  onto  the  original  finite  element  basis 
functions.  The  local  least  squares  projections,  however,  are  done  node  by  node  and  only 
consider  gradient  information  from  subdomains  of  Cl  which  contain  the  current  node.  Thus, 
a  series  of  smaller  projections  are  done  in  combination  with  some  averaging  techniques  to 
obtain  a  continuous  gradient  projection  for  the  entire  domain,  Cl.  The  details  of  each  of 
these  projection  techniques  will  be  described  further  in  Chapter  4. 

The  || Pex  ~V^\\  is  similarly  approximated  except  that  we  construct  a  continuous,  quadratic 
approximation  of  pex  using  local  projections  of  the  discontinuous  linear  finite  element  ap¬ 
proximations.  Once  this  is  done  the  L 2  norm  can  be  calculated  as  usual. 

We  now  return  to  the  issue  of  adaptive  gridding  to  briefly  describe  the  remeshing  strategy. 


Remeshing  Strategy 


Once  error  estimates  are  obtained  for  each  element,  say  e,,  a  new  mesh  density  (or  element 
size),  d,  is  calculated  which  requires  equidistribution  of  the  element  errors  across  all  the 
elements.  For  example  if  we  wish  to  reduce  the  error  in  each  element  by  a  factor  of  7,  then 
the  target  error,  ex,  for  an  element  in  the  new  mesh  can  be  given  by 


eT  = 


'ye 

Vn 


(3.109) 


r 
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where  e  is  the  error  over  the  entire  domain  and  N  is  the  number  of  elements.  If  one  assumes 
that  the  finite  element  method  is  of  order  k,  then  it  is  reasonable  to  write 

=  chk  (3.110) 

eT  =  cdk ,  (3.111) 

where  h  is  the  current  element  size  and  d  is  the  “ideal”  element  size  we  seek.  Clearly,  we 
have  assumed  that  we  are  in  the  asymptotic  range  of  the  finite  element  method  and  that 
the  convergence  constant,  c,  is  the  same  for  both  meshes.  This  may  or  may  not  be  the 
case.  Nevertheless,  the  system  (3.110)  -  (3.111)  can  be  solved  for  the  new  mesh  density,  d. 
obtaining 

/  7e  \1/k 

d  =  \«7w)  *'  (3112) 

The  new  element  size  computation  is  done  for  each  of  the  dependent  variables  in  the 
problem,  e.g.  velocity,  pressure,  and  an  ideal  mesh  density  is  obtained  for  each.  The  one 
used  for  the  generation  of  the  new  mesh  is  the  minimum  of  each  of  these. 

The  details  of  the  actual  redefinition  of  the  mesh  are  omitted  here.  We  refer  the  interested 
reader  to  [24]  and  [25]. 

Despite  the  rather  major  assumptions  made  above,  this  adaptive  remeshing  strategy  works 
remarkably  well.  The  strategy  has  been  verified  in  a  series  of  numerical  experiments  for  a 
variety  of  flow  types  (see  [34],  [33],  [24],  [29],  [25]).  We  use  this  adaptive  strategy  in  the 
problems  investigated  below. 


3.9  Some  Numerical  Results 

In  this  section,  we  apply  the  computational  techniques  outlined  above  to  the  two  problems 
presented  in  §3.7  above  and  investigate  convergence  of  the  approximate  solution  as  the 
mesh  is  refined. 


3.9.1  Flow  around  a  Cylinder 

Consider  first  the  cylinder  problem  presented  in  §3.7.1  .  We  again  assume  parabolic  inflow 
into  the  channel  as  follows: 

'  (l+q)(l-y2)  ' 

0 


u(-2  ,y;  q)  = 


-l<y<l. 


(3.113) 


Dawn  L.  Stewart 


Chapter  3.  2-D  Flow  Problems 


50 


The  outflow  boundary  condition  is  modified,  however,  since  applying  the  Dirichlet  boundary 
condition  at  the  outflow  causes  numerical  instability  (see  [17]).  A  free  boundary  condition 
is  used  which  requires 

~P(L,  y\  ?)n  +  r(u(L,  y;  q))  •  n  =  0,  1 
v(L,  y;q)=  0  / 

where  n  is  the  outward  normal  at  the  end  of  the  channel, 
are  applied  on  the  top,  bottom,  and  cylinder  sidewalls. 

As  noted  earlier,  we  wish  to  approximate  the  sensitivity  of  the  solution  with  respect  to 
the  inflow  parameter  q.  The  sensitivity  boundary  conditions  at  the  inflow  are  obtained  as 
before,  with 


-  1  <  y  <  1  (3.114) 

No-flow  and  no-slip  conditions 


s(~2,y;q)  = 


(1  -  y2) 

0 


-i<y<i. 


(3.115) 


The  computational  outflow  boundary  conditions  for  sensitivity  become 


-r(L,y,q)h  +  r(s(L,y;q))-n=  0, 
s  u(L,y,q)=  0 


} 


- 1  <  y  <  l. 


(3.116) 


Numerical  results  for  this  flow  problem  were  generated  using  the  approximation  techniques 
outlined  in  §3.8.  The  Reynolds  number  for  the  calculations  was  Re  =  100,  the  length  of 
the  channel  was  8  (L  =  6),  and  the  sensitivity  parameter  was  q  =  0.  Contour  plots  of 
the  u,  ^-velocity  fields  as  well  as  the  u,  ^-velocity  sensitivities  are  given  in  Figures  (3.6)  - 
(3.9).  In  Figure  (3.5),  the  initial  and  adapted  grids  are  shown.  It  is  clear  that  the  mesh 
is  refined  around  the  cylinder  and  in  areas  of  large  velocity  gradients.  This  gives  improved 
approximations  of  the  velocity  field  as  can  be  seen  in  Figures  (3.6)  and  (3.7).  Since  the 
mesh  refines  on  the  velocity  field,  it  is  convenient  that,  for  this  problem,  the  sensitivity 
flow  field  is  similar  to  the  velocity  field,  so  that  as  the  mesh  refines  we  obtain  improved 
sensitivity  approximations  as  well  (see  Figures  (3.8)-(3.9)).  It  is  important  to  note  that 
this  is  not  always  the  case  as  is  shown  in  [8].  The  code  has  been  modified  by  Jeff  Borggaard 
to  adapt  on  the  sensitivity  field  as  well  and  some  results  for  this  will  be  shown  later. 
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a)  Initial  Mesh  (1219  nodes/563  elements) 


b)  Second  Adapted  Mesh  (1439  nodes/677  elements) 


c)  Fourth  Adapted  Mesh  (19600  nodes/10818  elements) 


Figure  3.5:  Initial  and  Adapted  Meshes  for  Cylinder  Problem 
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Figure  3.6:  u- Velocity  Contours  for  Flow  around  a  Cylinder 
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b)  Second  Adapted  Mesh 


c)  Fourth  Adapted  Mesh 


-1-0.83-0.67  -0.5-0.33-0.17  0  0.17  0.33  0.5  0.67  0.83 

Figure  3.7:  v- Velocity  Contours  for  Flow  around  a  Cylinder 
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b)  Second  Adapted  Mesh 


c)  Fourth  Adapted  Mesh 


-1-0.83-0.67  -0 . 5-0^33-0 !  17  . . 0  0.17  0.33  0.5 ’"6T67  0.83  l 

Figure  3.9:  v- Velocity  Sensitivity  Contours  for  Flow  around  a  Cylinder 
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a)  Initial  Mesh  (625  Nodes/276  Elements) 


b)  Second  Adapted  Mesh  (1094  Nodes/505  Elements) 


c)  Fourth  Adapted  Mesh  (3962  Nodes/1901  Elements) 


Figure  3.10:  Initial  and  Adapted  Meshes  for  Bump  Problem 


3.9.2  Flow  over  a  Bump 


We  next  consider  the  flow  over  a  bump  problem  presented  in  §3.7.2.  The  outflow  boundary 
condition  for  the  computations  is  again  taken  to  be  a  free  boundary  condition  as  follows: 


~P(L,  V,  g)n  +  r(u(L,  y,  q))  -  n 
v(L,  y,  q) 


:?} 


0  <  y  <  3. 


(3.117) 


The  Reynolds  number  for  the  calculations  was  Re  =  100  and  the  length  of  the  channel 
was  L  —  8.  The  initial  and  adapted  grids  are  shown  in  Figure  (3.10).  Contour  plots  for 
u,  v  velocities  and  sensitivities  are  displayed  in  Figures  (3.11)  -  (3.12).  Note  that  the  mesh 
refines  in  the  area  of  the  bump  and  the  elements  are  allowed  to  become  larger  downstream 
in  the  channel  where  the  flow  is  again  quadratic.  The  bump  problem  is  similar  to  the 
cylinder  problem  in  that  the  sensitivities  are  largest  in  the  same  areas  where  the  mesh 
refines.  Thus  we  obtain  improved  approximations  for  sensitivities  as  we  refine  the  mesh  for 
the  flow. 
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So  far  we  have  developed  sensitivity  equations  for  both  1-D  and  2-D  parameter  dependent 
differential  equations.  We  have  obtained  good  numerical  approximations  for  the  continuous 
sensitivity  equations  using  an  adaptive  finite  element  technique.  In  the  1-D  problem,  there 
were  parameter  ranges  for  which  it  became  difficult  to  obtain  good  sensitivity  approxima¬ 
tions.  For  the  2-D  flow  problems,  we  observed  that  significant  refinement  was  needed  to 
obtain  “good”  sensitivity  approximations.  In  the  next  chapter,  we  look  at  a  projection 
technique  for  obtaining  better  derivative  approximations  for  the  state.  We  will  show  that 
these  derivative  approximations  improve  the  accuracy  of  the  sensitivity  calculations  and, 
in  addition,  stabilize  those  calculations,  over  larger  parameter  ranges. 
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a)  u-Velocity  Contours  -  Initial  Mesh 


b)  u-Velocity  Contours  -  Second  Adapted  Mesh 


c)  u-Velocity  Contours  -  Fourth  Adapted  Mesh 


d)  v-Velocity  Contours  -  Initial  Mesh 


e)  v-Velocity  Contours  -  Second  Adapted  Mesh 


f)  v-Velocity  Contours  -  Fourth  Adapted  Mesh 

-0.07  -0.054  -0.038  -0.022  -0.006  0.01  0.026  0.042  0.058  0.074  0.09 


Figure  3.11:  u, v-Velocity  Contours  for  Flow  over  a  Bump 
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a)  u-Velocity  Sensitivity  Contours  -  Initial  Mesh 


b)  u-Velocity  Sensitivity  Contours  -  Second  Adapted  Mesh 


c)  u-Velocity  Sensitivity  Contours  -  Fourth  Adapted  Mesh 


f)  v-Velocity  Sensitivity  Contours  -  Fourth  Adapted  Mesh 


-1.1  -0.93  -0.75  -0.58  -0.4  -0.23  -0.056  0.12  0.29  0.47  0.64 


Figure  3.12:  u,v- Velocity  Sensitivity  Contours  for  Flow  over  a  Bump 


Chapter  4 


Gradient  Approximations 


4.1  Improving  the  Gradient  Approximations 


The  goal  of  this  chapter  is  to  describe  how  projection  techniques  can  be  used  to  obtain 
better  sensitivity  approximations  by  obtaining  better  state  gradient  approximations.  There 
are  a  number  of  ways  to  do  this.  In  the  following  chapter,  we  analyze  a  global  and  local 
projection  technique  for  calculating  continuous  gradient  approximations  and  evaluate  their 
impact  on  providing  improved  sensitivity  values. 


4.2  1-D  Model  Problem 


Recall  that  the  approximate  sensitivity  equation  (2.30)  for  the  1-D  model  problem  has  the 
form 


60 
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Here,  gN(£)  is  a  piecewise  step  function  providing  an  approximation  of  the  spatial  deriva¬ 
tives  needed  in  (4.1).  Even  in  higher  dimensions  or  when  using  higher  order  basis  functions, 
the  finite  element  derivatives  will  be  discontinuous  across  element  faces.  As  in  the  error 
estimation  technique  described  above,  we  seek  to  replace  gN(£)  with  a  continuous  function 
which  we  will  denote  gN(£)-  We  will  now  describe  global  and  local  projection  schemes  for 
constructing  this  function. 


4.2.1  A  Global  Projection  Scheme 


In  its  simplest  form,  this  approach  replaces  the  discontinuous  piecewise  constant  function 
9N(0  by  its  projection  onto  the  space  of  piecewise  linear  splines  on  the  mesh  defined  by 
the  nodes  Thus,  we  consider  the  space 

5"  ={«(?):  S©  =  |>f«)},  (4.3) 

where  hf  (-)  are  the  hat  functions  defined  in  §2.5.  Observe  that  Stf  C  SN  and  SN  contains 
functions  with  non-zero  trace  on  the  boundary  of  [0, 1]. 


Here,  we  let  gN(€)  be  the  orthogonal  projection  of  gN(£)  onto  SN.  In  particular, 

n+i 

«"«)  =  £  ftAf  (f)  (4.4) 

i=0 

where  j3  =  (/?o,  A,  •  • .  ,  (3n+\)t  is  the  solution  of  a  linear  system  of  the  form 

MG3  =  Fa  (4.5) 

and  a  =  (al5  a2, . . .  ,  aN)T  contains  the  coefficients  defined  by  the  finite  element  approxi¬ 
mation  of  zN  (£).  Therefore, 


N 


9N( f)=Eai^Af(f), 

Mg  is  the  ( N  +  2)  x  ( N  +  2)  global  mass  matrix 

MG=[(hf(-),hf(-))]r, 

and  F  is  the  (N  +  2)  x  N  matrix 


F  = 


"(•).«•)' 1T 


for  i  —  1, 2, . . .  ,N  and  j,  k  =  0, 1, 2, . . .  ,  JV  +  1. 


(4.6) 


(4.7) 

(4.8) 
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4.2.2  A  Local  Projection  Scheme 

In  addition  to  the  global  projection  scheme,  we  consider  a  local  projection  scheme  which 
involves  performing  a  series  of  local  projections  on  subdomains  of  D  =  [0, 1],  At  each 
element  vertex,  we  define  the  subdomain  Qz  to  be  the  union  of  all  elements  for  which  & 
is  a  vertex.  In  define  </,(£)  to  be  the  least  squares  projection  of  ^(Oln  onto  the  space 
of  linear  polynomials  spanned  by  monomial  basis  functions.  Denote  these  basis  functions 
as  Pi{0  =  1  and  P2(£)  =  C  On  the  subdomain  S7j,  we  express  the  projection  as 

2 

9i(0  =  YlaiPi&'  (4.9) 

i—l 

where  the  vector  a  =  (ai,a2)T  contains  the  coefficients  of  the  basis  functions.  These 
coefficients  are  determined  by  solving  the  normal  equation  system 

MLa  =  b.  (4.10) 

The  matrix  Ml  is  of  the  form 

^=[<W,^(*)>lnJ>  (4-11) 

for  i,j  =  1,2.  The  vector  on  the  right  side  of  the  equation  is 

*=[<ao.sw(-)>U].  (4.i2) 

for  i  =  1,2.  Then,  on  Q,  we  define  the  continuous  local  projection  to  be 

_  JV+l 

9N(0  =  (4.13) 

t=0 

With  higher  order  finite  elements  and  in  higher  dimensions,  one  must  resolve  the  value  of 
9N(0  at  a  non-vertex  node.  An  averaging  technique  is  generally  used.  The  technique  will 
be  described  in  more  detail  in  the  next  section. 


4.2.3  Numerical  Results 

Derivative  Approximations 

Figures  (4.1)  and  (4.2)  show  the  finite  element  derivative  and  its  local  and  global  projections 
against  the  exact  spatial  derivative  for  a  N  =  4,  q  —  2.0  and  N  =  8,  q  =  1.2, 
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Derivative  and  Approximations  with  N  =  4  and  q  =  2 


Figure  4.1:  Finite  Element  Derivatives  with  Projections  at  N  =  4  and  q  =  2 


Derivative  and  Approximations  with  N  =  8  and  q  =  1.2 


Figure  4.2:  Finite  Element  Derivatives  with  Projections  at  N  =  8  and  q  =  1.2 
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L 2  Errors 

N  =  4  and  q  —  2 

N  =  8  and  q  =  1.2 

PWC 

0.4330 

1.1432 

Global 

0.2165 

0.8760 

Local 

0.2954 

1.0501 

Table  4.1:  L2  Errors  for  the  Derivative  Approximations 


respectively.  Clearly,  the  global  and  local  projections  give  different  results.  Since  we  have 
an  exact  solution,  and  thus  its  spatial  derivative,  we  can  calculate  element  by  element  L2 
errors  as  well  as  overall  L2  errors  for  the  whole  domain  for  each  derivative  approximation. 
Figures  (4.3)  -  (4.4)  show  the  L2  element  errors  for  the  two  cases,  q  —  2.0,  /N  =  4  and 
q  =  1.2,  /N  =  8,  repectively.  The  L 2  errors  over  the  whole  domain  for  each  of  these 
cases  is  summarized  in  Table  (4.1).  Note  that  the  errors  for  the  local  projection  technique 
are  highest  near  the  boundary,  but  that  away  from  the  boundary  the  local  projection 
technique  actually  has  less  error  than  the  global  projection  technique  for  these.  We  now 
analyze  these  techniques  in  calculating  sensitivities  for  varying  discretizations  (N,  M)  and 
parameter  values  {q)  to  understand  better  how  the  improved  derivative  approximations 
affect  our  numerical  sensitivities. 


Sensitivity  Approximations 


Figures  (4.5),  (4.6),  and  (4.7)  display  the  sensitivity  approximations  using  the  three  dif¬ 
ferent  derivative  approximations  against  the  exact  sensitivities  for  q  =  2.0,  1.4,  and  1.2, 
respectively.  It  is  clear  that  the  use  of  a  projection  technique  greatly  improves  the  accuracy 
of  our  sensitivity  approximations,  especially  as  q  ->•  1.  Recall  that  at  q  =  1.4  and  q  =  1.2 
our  sensitivity  approximations  obtained  using  PWC  derivatives  were  extremely  inaccurate 
so  that  the  approximations  do  not  even  show  up  on  the  graphs  in  Figures  (4.6)  and  (4.7) 
for  some  values  of  N.  As  shown  in  Figure  (4.8),  the  local  projections  clearly  decrease  the  L 2 
error  of  the  sensitivities  for  a  given  N  and  q.  Moreover,  the  most  promising  result  is  that 
the  projections  stabilize  the  calculations  over  the  parameter  range.  The  local  projections 
give  slightly  better  sensitivity  approximations  for  some  values  of  q.  Figure  (4.9)  compares 
the  L2  error  of  sensitivities  calculated  using  local  projections  with  the  error  calculated  using 
global  projections.  Note  that  as  N  increases,  the  global  projections  do  a  better  job  than 
local  projections  over  a  wider  range  of  q.  This  is  not  surprising  since  the  global  projec¬ 
tion  is  the  “best”  least  squares  linear  approximation  (using  the  finite  element  basis)  of  the 
piecewise  continuous  finite  element  derivative  in  the  limit. 
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Jig 


N=M=8  N=M=16 


Figure  4.5:  Sensitivity  Approximations  at  q  =  2 


N=M=4  N=M=8 


N=M=16 


N=W=32 


Figure  4.6:  Sensitivity  Approximations  at  q  =  1.4 
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N=M=8  N=M=16 


Figure  4.7:  Sensitivity  Approximations  at  q  =  1.2 


I  9  Prrnr  nf  .C^anertii/rti^c- 


q 


Figure  4.8:  Model  Problem  -  L2  Error  of  Sensitivity  Approximations  (PWC  and  Local) 
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Figure  4.9:  Model  Problem  -  L2  Error  of  Sensitivity  Approximations  (Global  and  Local) 
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Optimization  Results 

We  evaluated  the  same  two  cases  considered  in  §2.6.2  above:  (Case  1)  qinit  =  1.2  with 
q*  ~  2  and  (Case  2)  qinit  =  2.0  with  q*  ~  1.4.  This  time  the  simulations  were  performed 
using  the  piecewise  linear  derivative  approximations.  Tables  (4.2)  -  (4.5)  show  the  results 
of  these  simulations  a sN,M  ranged  from  2- 128.  Note  that  the  use  of  the  piecewise  linear 
derivative  approximations  clearly  improved  the  results  of  the  optimization  algorithm  for 
Case  1  and  that  the  global  and  local  projection  schemes  provided  very  similar  results.  The 
improvement  was  not  as  marked  for  Case  2,  however  the  scheme  did  converge  for  the  case 
N  —  M  =  16  with  the  improved  gradient  approximations  and  did  not  with  the  piecewise 
constant  derivatives.  In  addition,  the  local  scheme  to  slightly  less  time  to  converge  than 
the  global  projection  scheme  (see  Tables  (4.4)  -  (4.5)). 
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N 

M 

CONV/DNC 

ITER 

TIME 

J"|| 

q 

2 

2 

CONV 

20 

9.43 

.32586 

1.9604 

2 

5 

CONV 

16 

9.06 

.32548 

1.9514 

4 

4 

CONV 

15 

17.80 

.26463 

1.9509 

4 

9 

CONV 

15 

12.91 

.26468 

1.9481 

8 

8 

CONV 

12 

27.01 

.24327 

1.9423 

8 

17 

CONV 

13 

25.84 

.24326 

1.9421 

16 

16 

CONV 

13 

27.95 

.24235 

1.9375 

16 

33 

CONV 

13 

31.62 

.24235 

1.9376 

32 

32 

CONV 

12 

77.32 

.24272 

1.9361 

32 

65 

CONV 

12 

78.90 

.24273 

1.9361 

64 

64 

CONV 

12 

288.06 

.24274 

1.9359 

64 

129 

CONV 

11 

302.73 

.24274 

1.9359 

128 

128 

CONV 

11 

355.10 

.24285 

1.9358 

Table  4.2:  Optimization  Results  for  Case  1,  Global  Projection  Scheme 


N 

M 

CONV/DNC 

ITER 

TIME 

\\jn\\ 

q 

2 

2 

CONV 

20 

9.43 

.32586 

1.9604 

2 

5 

CONV 

16 

9.06 

.32548 

1.9514 

4 

4 

CONV 

15 

17.80 

.26463 

1.9509 

4 

9 

CONV 

15 

12.91 

.26468 

1.9481 

8 

8 

CONV 

12 

27.01 

.24327 

1.9423 

8 

17 

CONV 

13 

25.84 

.24326 

1.9421 

16 

16 

CONV 

13 

27.95 

.24235 

1.9375 

16 

1  33 

CONV 

13 

31.62 

.24235 

1.9376 

32 

32 

CONV 

12 

77.32 

.24272 

1.9361 

32 

65 

CONV 

12 

78.90 

.24273 

1.9361 

64 

64 

CONV 

12 

288.06 

.24274 

1.9359 

64 

129 

CONV 

11 

302.73 

.24274 

1.9359 

128 

128 

CONV 

11 

355.10 

.24285 

1.9358 

Table  4.3:  Optimization  Results  for  Case  1,  Local  Projection  Scheme 
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N 

M 

CONV/DNC 

ITER 

TIME 

q 

2 

2 

DNC 

20 

2.06 

.93423 

2.3054 

2 

5 

DNC 

20 

64.37 

.72262 

1.3530 

4 

4 

DNC 

20 

133.51 

.53248 

1.4557 

4 

9 

DNC 

20 

132.27 

.52703 

1.4405 

8 

8 

DNC 

20 

127.18 

.48242 

1.5029 

8 

17 

DNC 

20 

143.56 

.50490 

1.4791 

16 

16 

CONV 

14 

28.69 

.41967 

1.4329 

16 

33 

CONV 

15 

33.07 

.41973 

1.4322 

32 

32 

CONV 

!  15 

13.22 

.41001 

1.4197 

32 

65 

CONV 

11 

12.75 

.41001 

1.4196 

64 

64 

CONV 

13 

45.43 

.40819 

1.4168 

64 

129 

CONV 

14 

48.74 

.40819 

1.4167 

128 

128 

CONV 

13 

258.93 

.40774 

1.4160 

Table  4.4:  Optimization  Results  for  Case  2,  Global  Projection  Scheme 


N 

M 

CONV/DNC 

ITER 

TIME 

ik"  ii 

q 

2 

2 

DNC 

20 

60.82 

.73574 

1.4154 

2 

5 

DNC 

20 

64.60 

.72273 

1.3634 

4 

4 

DNC 

20 

130.65 

.53216 

1.4558 

4 

9 

DNC 

20 

131.56 

.52991 

1.4513 

8 

8 

DNC 

20 

133.27 

.50451 

1.4830 

8 

17 

DNC 

20 

133.92 

.50470 

1.4807 

16 

16 

CONV 

14 

30.00 

.41968 

1.4328 

16 

33 

CONV 

16 

32.03 

.41974 

1.4322 

32 

32 

CONV 

11 

11.21 

.41000 

1.4197 

32 

65 

CONV 

15 

14.34 

.41001 

1.4196 

64 

64 

CONV 

13 

44.45 

.40820 

1.4168 

64 

129 

CONV 

12 

45.70 

.40819 

1.4168 

128 

128 

CONV 

12 

253.86 

.40774 

1.4160 

Table  4.5:  Optimizations  Results  for  Case  2,  Local  Projection  Scheme 
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Figure  4.10:  Typical  Subdomain  of  an  Element  Vertex  £ 

4.3  A  Local  Projection  for  Higher  Dimensions 


In  this  section,  we  describe  the  local  projection  scheme  employed  in  the  error  estimation 
module  of  the  finite  element  code  described  in  Chapter  2.  This  local  projection  scheme  is 
virtually  the  same  as  the  local  scheme  for  the  1-D  Model  Problem  discussed  in  Section  4.2.2 
above,  yet  is  presented  for  a  two  dimensional  problem  below  since  the  use  of  higher  order 
finite  elements  in  higher  dimensions  adds  some  complexities  which  were  not  previously 
discussed. 

Again,  the  local  projection  scheme  involves  performing  a  series  of  projections  on  subdomains 
of  D.  The  projected  gradients  are  given  as  polynomial  expansions  around  a  given  vertex, 
£»?  of  a  finite  element  mesh.  The  subdomain,  over  which  the  projection  consists  of 
all  elements  having  £t-  =  (xi:  y{)  as  a  vertex.  Figure  (4.10)  illustrates  a  typical  subdomain 
a  finite  element  mesh  of  quadratic  triangles.  In  principle,  the  choice  of  the  degree  of 
the  polynomial  expansion  for  the  improved  gradient  approximation  is  independent  of  the 
selection  of  the  finite  element  basis  being  used.  However,  in  practice,  the  degree  of  the 
polynomial  expansion  is  chosen  to  match  the  degree  of  the  finite  element  basis  employed. 
This  leads  to  an  order  of  accuracy  improvement  in  the  gradient  approximations.  For  all 
the  numerical  results  presented  in  Chapter  5  below,  quadratic  triangular  finite  elements 
were  used  and  so  the  locally  projected  gradients  are  written  as  polynomials  of  degree  two 
on  . 


Dawn  L.  Stewart 


Chapter  4.  Gradient  Approximations 


73 


At  element  vertices,  ^2,  and  £3,  we  define  g*  to  be  the  local  least  squares  projection  of 
the  finite  element  derivative,  Vuh,  onto  the  space  of  quadratic  polynomials.  For  ease  of 
notation,  we  denote  VuA  by  gk.  Letting  P  =  [1,  x,  y,  x2,  x y,  y2}  denote  the  basis  functions 
of  this  space,  we  can  express  each  component  of  the  gradient  projection,  g*  and  g*,  as 

g*  =  P’a*  and,  g*  =  Praj,  (4.14) 

where  the  vectors  a*,  By  €  1R6  contain  the  coefficients  of  the  basis  functions.  These  coeffi¬ 
cients  are  obtained  by  solving  the  following  least-squares  problems: 

mir4  /  fe  -  Sx)  and,  min^  [  (g*  -  gj)  dQ  (4.15) 

Jau  z 

Thus,  for  each  component  of  g*,  we  solve  the  following  6x6  system  of  linear  equations  for 
and  ay . 


pTPdn 


{a*}  = 


PTP  dQ 


K}  = 


Prgirfn| 

(4.16) 

PTg£<*«  |  ■ 

(4.17) 

The  finite  element  fluxes,  g£  and  g £,  are  obtained  in  the  usual  manner  by  differentiating 
the  finite  element  basis  functions.  Note  that  the  left  hand  matrix  is  independent  of  the 
quantity  being  projected  and  thus  can  be  viewed  as  the  projection  matrix  for  node,  £f, 
and  can  be  used  for  obtaining  locally  projected  derivative  approximations  for  all  of  the 
dependent  variables  (e.g.  u,  v,  and  T)  as  long  as  the  projection  basis,  P,  is  not  changed. 


Once  we  are  done,  we  have  a  quadratic  expression  for  the  locally  projected  derivative,  g*, 
at  each  vertex.  Let  g^  denote  the  expression  for  g*  obtained  by  solving  the  systems  (4.16)- 
(4.17)  where  is  the  subdomain  associated  with  element  vertex,  Define  g|2  and  g£3 
similarly  for  the  remaining  element  vertices,  £2  and  £3.  We  need  a  unique  definition  of  g* 
for  any  point,  £  inside  the  element  (see  Figure  (4.11)).  For  quadratic  elements,  there  are 
several  ways  to  do  this.  I  describe  the  technique  employed  in  the  current  version  of  the 
code. 


Unique  nodal  values  of  g*  at  the  element  vertices,  which  we  denote  gj,  g£,  and  g3,  are 
simply  defined  as  follows: 

«?  =  «&(&),  *  =  1,2,3.  (4.18) 

Nodal  values  for  the  midside  nodes  are  obtained  by  averaging  the  values  of  the  polynomial 
expressions  for  g*  at  the  endpoints  of  element  side.  For  example, 

§*(£12)  =  2  (g|i(£l2)  +g|2(£l2))  • 


(4.19) 
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Figure  4.11:  Element  with  Three  Quadratic  Expressions  for  g* 

Then,  at  any  point  £  in  the  element,  the  value  of  the  locally  projected  derivative  g*  at  £  is 

6 

g  *(«)  =  I>;g*  (4.20) 

j= 1 

where  the  Nj  are  the  quadratic  basis  functions  for  the  finite  element  space  and  the  g*  are 
the  nodal  values  of  the  locally  projected  derivative  obtained  as  describe  above. 

In  the  next  chapter,  we  apply  the  local  projection  technique  described  above  in  order  to 
obtain  improved  sensitivity  approximations  for  two  flow  problems. 


Chapter  5 

Numerical  Results  for  2-D  Problems 


We  now  return  to  the  two  specific  flow  problems  described  in  §3.7  to  see  if  the  projection 
techniques  described  in  §4.3  can  be  used  to  obtain  improved  sensitivity  approximations. 
We  begin  by  considering  the  cylinder  problem  discussed  in  §3.7.1. 


5.1  Flow  Around  a  Cylinder 


Numerical  approximations  to  the  state  and  sensitivities  for  this  problem  were  calculated 
over  a  range  of  Reynolds  numbers.  In  each  case,  the  results  obtained  for  the  initial  and 
first  adapted  meshes  were  compared  to  an  approximation  which  was  generated  by  adapting 
to  a  very  fine,  final  mesh  and  then  interpolating  the  “true”  solution  from  the  final  mesh 
onto  the  initial  or  first  adapted  meshes.  On  this  final  mesh,  both  schemes  converged  to 
solutions  differing  by  less  than  10-3.  For  comparison  purposes,  we  chose  the  unprojected 
solution  on  the  final  mesh  as  the  “true”  solution. 

We  note  that  the  length  of  the  channel  was  8  (i.e. ,  L  =  6)  for  all  of  the  runs  in  this  section. 
This  length  was  not  sufficient  for  the  flow  at  the  outflow  to  return  to  the  parabolic  inflow 
especially  over  the  range  of  Reynolds  numbers  being  considered,  however,  it  was  sufficient 
to  meet  the  computational  “free”  outflow  boundary  condition.  We  will  use  the  Re  =  350 
case  to  show  that  the  results  for  the  higher  Reynolds  cases,  were  not  affected  by  the  length 
of  the  channel. 


In  this  chapter,  we  will  also  use  an  adaptive  technique  which  refines  on  sensitivity  errors  as 
well  as  flow  errors.  This  technique  is  completely  analogous  to  the  technique  used  to  refine 
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on  the  flow  errors  (see  [9]).  In  the  results  which  follow,  we  will  be  careful  to  identify  meshes 
which  were  adapted  on  approximations  of  flow  errors  and  meshes  which  were  adapted  on 
approximations  of  both  flow  and  sensitivity  errors. 

We  present  a  detailed  error  analysis  of  the  velocity  sensitivities  s  for  two  Reynolds  numbers, 
Re  =  100  and  Re  —  350.  We  begin  with  the  results  at  Re  =  100.  Figure  (5.1)  shows  the 
initial  mesh,  the  first  adapted  meshes,  and  final  mesh  for  this  case.  Notice  that  we  present 
three  first  adapted  meshes:  one  adapted  on  the  flow  only,  one  adapted  on  both  the  flow 
and  the  sensitivities  calculated  with  unprojected  derivatives,  and  one  adapted  on  both  the 
flow  and  sensitivities  calculated  with  projected  derivatives.  At  this  Reynolds  number,  the 
difference  between  these  three  meshes  is  less  dramatic  than  it  is  at  Re  =  350,  yet  we  present 
all  of  them  for  consistency  and  comparison  purposes. 

As  expected,  the  scheme  employing  the  locally  projected  derivatives  gives  better  sensitivity 
results  as  seen  in  Figures  (5.2)  and  (5.3).  In  these  figures,  recall  that  the  “true”  solution 
is  the  unprojected  solution  on  the  final  mesh  interpolated  onto  the  initial  mesh.  Once  this 
interpolation  is  made,  node  by  node  errors  can  be  calculated.  These  errors  are  shown  in 
Figure  (5.4).  From  these  plots,  we  can  see  that  the  local  projection  scheme  reduces  the 
maximum  error  by  about  a  factor  of  2. 

We  also  look  at  the  sensitivity  solutions  using  schemes  on  the  first  adapted  mesh.  The 
first  adapted  mesh  (for  the  flow)  is  shown  in  Figure  (5.1b).  The  sensitivity  results  for  the 
unprojected  and  locally  projected  derivative  schemes  for  this  mesh  are  shown  in  Figures 
(5.5)  and  (5.6).  In  these  figures,  the  “true”  solution  is  now  the  unprojected  solution  on  the 
final  mesh  interpolated  onto  the  first  adapted  mesh.  Note  that  the  difference  between  the 
approximations  using  the  projected  gradients  and  those  using  the  “natural”  finite  element 
derivatives  is  no  longer  as  apparent.  Again,  a  node  by  node  comparison  was  made.  These 
errors  are  shown  in  Figures  (5.7a, b)  -  (5.8a, b)  and  were  plotted  on  the  same  scales  used 
for  the  initial  mesh  for  easy  comparison.  We  also  plotted  nodal  errors  for  the  two  schemes 
in  the  cases  where  we  adapted  on  the  flow  and  the  sensitivities  (see  Figures  (5.7c,d)  - 
(5.8c, d)).  At  this  Reynolds  number,  there  does  not  seem  to  be  a  significant  advantage 
gained  by  adapting  on  both  the  flow  and  sensitivity  errors. 
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a)  Initial  Mesh  (1219  nodes/563  elements) 


b)  First  Adapted  Mesh  -  Flow  Only  (956  nodes/444  elements) 


c)  First  Adapted  Mesh  -  Flow  &  Unprojected  Sens  (1 082  nodes/504  elements) 


d)  First  Adapted  Mesh  -  Flow  &  Projected  Sens  (1 066  nodes/496  elements) 


Figure  5.1:  Meshes  for  Cylinder  Problem  at  Re  =  100 
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a)  Unprojected 


b)  Projected 


-1  -0.83  -0.67  -0.5  -0.33  -0.X7  0  0.17  0.33  0.5  0.67  0.83  1 

Figure  5.3:  v- Velocity  Sensitivities  on  Initial  Mesh  for  Re  =  100 
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a)  Error  in  u  Sensitivitities  -  Unprojected  Derivatives 

0  0.082  0.16  0.24  0.33  0.41  0.4$ 


b)  Error  in  u  Sensitivitities  -  Projected  Derivatives 


c)  Error  in  v  Sensitivitities  -  Unprojected  Derivatives 


a^mmmmmmmmmmmFKrrrTrrrrr^r"-  v  '  v .  ■  ■f^.Kn^'irsaaaaaaaaaK 

0  0.085  0.17  0.25  0.34  0.42  0.51 


d)  Error  in  v  Sensitivitities  -  Projected  Derivatives 


Figure  5.4:  Error  of  Sensitivity  Approximations  on  Initial  Mesh  for  Re 
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b)  Projected  Derivatives 


ascatr  TfcsKTs rr 
0.092 


Figure  5.5:  u- Velocity  Sensitivities  on  First  Adapted  Mesh  for  Re  =  100 
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a)  Unprojected  Derivatives 


b)  Projected  Derivatives 


r  ™ 1 . :  ~rzr.  ^  ■  : - : . ■  .■ 


-1  -0.83  -0.67  -0.5  -0.33  -0.17  0  0.17  0.33  0.5  0.67  0.83  1 


Figure  5.6:  v- Velocity  Sensitivities  on  First  Adapted  Mesh  for  Re  —  100 
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c)  Unprojected  Derivatives,  Adapt  on  Flow  and  Sensitivity 


0  0.082  0.16  0.24  0.33  0.41  0.49 


Figure  5.7:  Error  of  u- Velocity  Sensitivity  Approximations  on  First  Adapted  Mesh  for 
Re  =  100 


Dawn  L.  Stewart 


Chapter  5.  Numerical  Results  for  2-D  Problems 


84 


c)  Unprojected  Derivatives,  Adapt  on  Flow  and  Sensitivity 


0 


0.34 


0.42  0.51 


Figure  5.8:  Error  of  v- Velocity  Sensitivity  Approximations  on  First  Adapted  Mesh  for 
Re  =  100 
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Initial  Mesh 

Unproj 

Proj 

u 

V 

Su 

sv 

1.7397E-01 

9.0160E-02 

5.0463E-01 

2.2453E-01 

1.7397E-01 

9.0160E-02 

2.8227E-01 

1.3050E-01 

First  Adapted  Meshes 

Unproj  (Flow) 

Proj  (Flow) 

Unproj  (Flow  Sc  Sens) 

Proj  (Flow  Sc  Sens) 

u 

1.0313E-01 

1.0313E-01 

7.2637E-02 

5.3305E-02 

V 

4.5997E-02 

4.5997E-02 

3.3000E-02 

2.8115E-02 

3.4173E-01 

3.5800E-01 

2.5988E-01 

2.7273E-01 

Sv 

1.4487E-01 

1.1039E-01 

9.1075E-02 

7.6292E-02 

Table  5.1:  L2  Errors  for  Flow  and  Sensitivities  at  Re  =  100 


In  order,  to  get  an  overall  evaluation  of  the  errors,  we  used  the  node  by  node  errors  to 
calculate  an  L2  error  over  the  entire  domain  for  the  flow  and  the  sensitivities  (see  Table 
(5.1)).  This  more  clearly  shows  the  50%  reduction  of  the  error  for  the  local  projection 
scheme  on  the  initial  mesh.  Note  that  the  errors  for  u,  v  on  the  first  adapted  meshes  are 
the  same  for  the  unprojected  and  projected  schemes  on  the  mesh  that  was  refined  only  on 
the  flow,  but  are  slightly  different  (since  the  meshes  are  slightly  different)  when  the  meshes 
were  also  refined  on  the  sensitivity.  Note  that  adapting  on  the  flow  and  sensitivity  errors 
not  only  improved  the  sensitivity  approximations  over  those  obtained  by  adapting  on  the 
flow  alone,  but  also  improved  the  numerical  approximations  for  the  flow. 

We  now  turn  to  the  results  for  Re  =  350.  We  begin  by  comparing  the  state  and  sensitivity 
approximations  for  L  =  6  and  L  =  15  to  ascertain  whether  or  not  the  length  of  the  channel 
affects  our  results.  Figure  (5.9)  shows  the  initial  meshes  for  the  two  channel  lengths  as  well 
as  the  u,  v  velocity  contours.  It  is  clear  that  the  flow  more  nearly  returns  to  the  inflow  for 
the  longer  channel.  However,  the  results  for  the  shorter  channel  are  very  similar  to  those  for 
the  longer  channel  in  the  areas  where  they  overlap.  In  fact,  Figures  (5.10)  and  (5.11)  show 
that  the  dramatic  difference  between  the  sensitivities  obtained  with  unprojected  derivatives 
and  those  obtained  using  locally  projected  derivatives  occur  in  both  the  long  and  the  shorter 
channels.  This  gives  us  reasonable  confidence  that  the  length  of  the  channel  is  not  affecting 
our  results. 

We  now  return  to  complete  the  same  error  analysis  for  Re  =  350  that  we  did  for  Re  = 
100.  The  initial  and  first  adapted  meshes  are  shown  in  Figure  (5.12).  Note  that  the 
difference  between  the  three  first  adapted  meshes  is  definitely  greater  at  this  Reynolds 
number.  This  is  due,  at  least  in  part,  to  the  greater  discrepancy  between  the  flow  error, 
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the  unprojected  sensitivity  error,  and  the  projected  sensitivity  error.  At  this  Reynolds 
number,  the  differences  between  the  sensitivities  calculated  with  the  two  different  derivative 
schemes  is  now  dramatic  as  can  be  seen  in  Figures  (5.13)  -  (5.14).  The  node  by  node  error 
analysis  (see  Figure  (5.15))  shows  that  locally  projected  derivatives  are  definitely  better, 
although  at  this  Reynolds  number,  both  approximations  contain  fairly  larger  errors.  Using 
the  locally  projected  derivatives  to  calculate  sensitivities  reduces  the  overall  L2  error  by 
about  600%  (see  Table  (5.2)). 

The  same  analysis  is  done  on  the  first  adapted  meshes  (See  Figures  (5.16)  -  (5.19)).  Here, 
adapting  on  the  flow  as  well  as  the  sensitivity  makes  a  greater  difference  in  the  error 
reduction  from  the  initial  mesh  to  the  first  adapted  mesh.  This  is  easily  seen  in  the 
overall  L2  errors  for  the  approximate  flow  and  sensitivities  displayed  in  Table  (5.2).  It  is 
interesting  to  note  that  the  errors  for  all  the  quantities  are  less  for  the  Projected  (Flow  & 
Sensitivity)  than  they  are  for  the  Unprojected  (Flow  &  Sensitivity)  despite  the  fact  that 
the  Unprojected  (Flow  &  Sensitivity)  mesh  is  quite  a  bit  finer. 

To  conclude,  we  observe  that  using  the  locally  projected  derivatives  clearly  stabilizes  the 
calculations  over  a  larger  range  of  Reynolds  numbers  as  it  did  in  §4.2.3  for  the  1-D  model 
problem.  Plots  of  the  overall  L2  sensitivity  errors  on  the  initial  mesh  are  shown  in  Figures 
(5.20)  -  (5.21).  In  addition,  the  mesh  refinement  very  effectively  reduces  the  errors  for 
the  flow  and  the  sensitivities.  At  most  Reynolds  numbers,  after  the  mesh  is  refined  once, 
there  is  little  difference  between  the  sensitivities  calculated  with  the  unprojected  and  lo¬ 
cally  projected  techniques.  At  higher  Reynolds  numbers,  however,  using  locally  projected 
derivatives  to  calculate  sensitivities  remains  advantageous. 

We  now  turn  our  attention  to  obtaining  numerical  approximations  to  the  design  sensitivities 
for  the  flow  over  a  bump  problem. 
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Initial  Mesh 

Unproj 

Proj 

u 

V 

Su 

Sty 

1.4808E-00 

4.1024E-01 

1.9438E+01 

3.0582E-00 

1.4808E-00 

4.1024E-01 

3.6791E-00 

5.8348E-01 

First  Adapted  Meshes 

Unproj  (Flow) 

Proj  (Flow) 

Unproj  (Flow  &  Sens) 

Proj  (Flow  &  Sens) 

u 

9.6580E-01 

9.6580E-01 

3.5985E-01 

2.8674E-01 

V 

2.6044E-01 

2.6044E-01 

1.0257E-01 

7.7747E-02 

3.8080E-00 

3.1614E-00 

1.1145E-00 

1.0423E-00 

Sty 

1.5339E-00 

8.7266E-01 

3.4686E-01 

2.9430E-01 

Table  5.2:  L2  Errors  for  Flow  and  Sensitivities  at  Re  =  350 
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b)  Initial  Mesh  -  L  =  15 


c)  u  Velocity  Contours  -  L  =  6 
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d)  u  Velocity  Contours  -  L  =  1 5 


V  7  u 

'U-r^%r 


BQ-- 


e)  v  Velocity  Contours  -  L  =  6 
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f)  v  Velocity  Contours  -  L  =  1 5 


Figure  5.9:  Initial  Meshes  and  u,v- Velocity  Contours  for  L  =  6, 15  and  Re  =  350 
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a)  u  Velocity  Sensitivity  Contours,  Unprojected  Derivatives:  L  =  6 

— liimil  i  "  ~im . . r~* - -  wvtnfHiiawi 
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b)  u  Velocity  Sensitivity  Contours,  Projected  Derivatives:  L  =  6 


c)  u  Velocity  Sensitivity  Contours,  Unprojected  Derivatives:  L  =  15 

-0.95  -0.53  -0.1  0.32  0.74  1.2  1.6 


d)  u  Velocity  Sensitivity  Contours,  Projected  Derivatives:  L  =  15 

Figure  5.10:  u-Velocity  Sensitivity  Contours  for  L  =  6,15  and  Re  =  350 
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a)  v  Velocity  Sensitivity  Contours,  Unprojected  Derivatives:  L  =  6 
-TSTm -(L 48  -0.32  -0.16  0  0.16  0 . 32  OTo^sT*' 0 . 8 


b)  v  Velocity  Sensitivity  Contours,  Projected  Derivatives:  L  =  6 


c)  v  Velocity  Sensitivity  Contours,  Unprojected  Derivatives:  L  =  15 


0.64  0.8 


d)  v  Velocity  Sensitivity  Contours,  Projected  Derivatives:  L  =  15 

Figure  5.11:  v- Velocity  Sensitivity  Contours  for  L  =  6, 15  and  Re  =  350 
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b)  First  Adapted  Mesh  -  Flow  Only  (1164  nodes/542  elements) 


c)  First  Adapted  Mesh  -  Flow  &Unprojected  Sens  (2249  nodes/1071  elements) 


d)  First  Adapted  Mesh  -  Flow  &  Projected  Sens  (1556  nodes/732  elements) 


Figure  5.12:  Meshes  for  Cylinder  Problem  at  Re  =  350 


c)  "True"  Solution 


-0.95  -0.53  ' ‘-671  0.32 

Figure  5.13:  u-Velocity  Sensitivities 
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Figure  5.16:  u- Velocity  Sensitivities  on  First  Adapted  Mesh  for  Re  =  350 
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a)  Unprojected  Derivatives 
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Figure  5.17:  v- Velocity  Sensitivities  on  First  Adapted  Mesh  for  Re  =  350 
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c)  Unprojected  Derivatives,  Adapt  on  Flow  and  Sensitivity 


d)  Projected  Derivatives,  Adapt  on  Flow  and  Sensitivity 


~ 0 " '  0.46  0.92  'l'l . :  1  ."8  2.3 .  '2.8 

Figure  5.18:  Error  of  u-Velocity  Sensitivity  Approximations  on  First  Adapted  Mesh  for 
Re  —  350 
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Figure  5.20:  Flow  About  Cylinder  -  L2  Error  of  u-Sensitivity  Approximations  on  Initial 
Mesh 


Figure  5.21:  Flow  About  Cylinder  -  L2  Error  of  v-Sensitivity  Approximations  on  Initial 
Mesh 
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5.2  Flow  over  a  Bump 


We  begin  by  doing  a  qualitative  comparison  of  the  sensitivity  values  we  obtain,  with  those 
presented  in  [11].  We  will  point  out  some  of  the  difficulties  of  calculating  shape  sensitivities 
and  show  how  the  process  of  error  estimation  and  grid  refinement  is  extremely  important 
in  obtaining  accurate  numerical  approximations  of  the  sensitivities  for  these  problems. 

Figure  (5.22)  displays  vector  sensitivity  plots  for  a  =  £*i  =  0.5,  A  =  0.5,  and  L  =  10.  These 
plots  are  qualitatively  very  comparable  to  those  presented  on  page  172  in  [11].  As  Burkardt 
notes,  the  shape  sensitivities  for  smaller  Reynolds  numbers  appear  as  “whirlpools”  and  are 
predominately  localized  to  the  region  above  the  bump.  As  the  Reynolds  number  increases, 
however,  the  effect  of  the  bump  is  carried  downstream.  This  is  especially  apparent  in  Figure 
(5.22c). 

Another  effect  of  increasing  the  Reynolds  number  on  flow  calculations,  is  that  the  task 
of  meeting  the  outflow  boundary  conditions  becomes  more  challenging  numerically.  We 
clearly  see  this  effect  in  the  Re  -  500  case  shown  in  Figure  (5.23).  With  a  channel  length 
°f  L  —  10,  the  flow  was  not  able  reach  the  parabolic  flow  profile  of  the  inflow.  So  for  this 
case,  we  lengthened  the  channel  to  L  =  20.  Notice  also  that  we  are  getting  a  large  secondary 
^whirlpool  further  downstream  from  the  bump.  This  phenomenon  did  not  appear  in  the 
sensitivities  presented  in  [11].  We  will  examine  this  further  for  the  case  Re  =  1000. 

We  used  the  case  of  Re  =  1000  to  examine  a  number  of  issues  relating  to  our  sensitivity 
approximations.  First,  we  investigate  the  secondary  “whirlpool”  which  appears  in  our 
sensitivity  calculations.  We  also  evaluate  the  accuracy  of  the  sensitivity  approximations 
near  the  bump.  We  do  this  by  using  the  error  estimation  and  grid  refinement  process 
presented  in  §3.8.4.  The  mesh  was  adapted  only  on  flow  errors.  Figure(5.24)  shows  the 
initial  mesh,  which  is  similar  in  density  to  the  one  used  in  [11],  and  the  adapted  meshes. 
Note  that  the  mesh  refines  where  one  would  expect,  in  the  regions  of  large  velocity  gradients 
around  and  downstream  from  the  bump.  It  is  also  refining  at  the  outflow  in  an  effort  to 
accurately  meet  the  outflow  boundary  condition. 

Figures  (5.25)  and  (5.26)  display  u  and  v  contours  for  the  flow  on  the  initial  and  adapted 
meshes.  It  is  important  to  note  that  as  the  mesh  refines,  the  contours  smooth  and  the 
accuracy  of  the  gradients  in  u  and  v  are  greatly  improved  especially  in  the  vicinity  of  the 
bump.  This  lack  of  accuracy  of  the  velocity  gradients  on  the  initial  mesh  greatly  hampers 
our  ability  to  obtain  good  sensitivity  approximations.  This  is  seen  by  noting  that  the 
sensitivity  boundary  conditions  on  the  bump,  (3.66)  -  (3.67)  require  approximations  to  the 
velocity  gradients  on  the  boundary.  If  there  are  large  errors  in  the  velocity  gradients,  there 
will  be  large  errors  in  the  sensitivity  approximations.  This  effect  can  be  seen  in  Figures 
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(5.27)  -  (5.29).  As  the  gradient  approximations  are  improved  in  the  second  and  third 
adapted  meshes,  the  values  of  the  sensitivities  become  much  more  accurate,  not  only  in  the 
local  area  of  the  bump,  but  downstream  as  well.  Note  also  that  the  secondary  “whirlpool” 
which  appeared  in  the  Re  =  500  case,  is  seen  on  the  initial  mesh  in  this  case  also.  As  the 
mesh  is  refined,  however,  the  size  of  the  “whirlpool”  is  diminished  until  it  is  almost  gone 
as  can  be  seen  in  Figure  (5.29c). 
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c)  Re  =  100 

Figure  5.22:  Sensitivity  Vectors,  L  =  10,  Flow  over  a  Bump 
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c)  u  Contour  Plot,  Re  =  500,  L  =  20 


d)  Sensitivity  Vector  Plot,  Re  =  500,  L  =  20  (Part  of  Channel) 


Figure  5.23:  u- Velocity  Contours  and  Sensitivity  Vectors,  Re  =  500,  Flow  over  a  Bump 
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a)  Initial  Mesh 


b)  Rrst  Adapted  Mesh 


c)  Second  Adapted  Mesh 


d)  Third  Adapted  Mesh 


Figure  5.24:  Flow  over  a  Bump,  Initial  and  Adapted  Meshes  for  Re  =  1000,  L  -  20 
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a)  u  Contour  Plot,  Re  =  1000,  Initial  Mesh 


b)  u  Contour  Plot,  Re  =  1000,  First  Adapted  Mesh 


c)  u  Contour  Plot,  Re  =  1 000,  Second  Adapted  Mesh 


d)  u  Contour  Plot,  Re  =  1000,  Third  Adapted  Mesh 
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Figure  5.25:  u-Velocity  Contours  for  Flow  over  a  Bump 
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b)  u  Sensitivity  Contour  Plot,  Re  =  1000,  First  Adapted  Mesh 


c)  u  Sensitivity  Contour  Plot,  Re  =  1000,  Second  Adapted  Mesh 


d)  u  Sensitivity  Contour  Plot,  Re  =  1000,  Third  Adapted  Mesh 
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Figure  5.27:  u- Velocity  Sensitivity  Contours  for  Flow  over  a  Bump 
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a)  v  Sensitivity  Contour  Plot,  Re  =  1 000,  Initial  Mesh 


d)  v  Sensitivity  Contour  Plot,  Re  =  1 000,  Third  Adapted  Mesh 
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Figure  5.28:  v- Velocity  Sensitivity  Contours  for  Flow  over  a  Bump 
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a)  Sensitivity  Vector  Plots,  Re  =  1000,  Initial  Mesh 


b)  Sensitivity  Vector  Plots,  Re  =  1000,  First  Adapted  Mesh 


c)  Sensitivity  Vector  Plots,  Re  =  1000,  Second  Adapted  Mesh 
Figure  5.29:  Sensitivity  Vector  Plots  on  Initial  and  Adapted  Meshes 
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Analysis  of  Cost  Functionals  and  Gradients  -  An  Optimization  Issue 

We  now  turn  to  the  issue  of  using  numerical  approximations  of  the  state  and  sensitivities 
to  approximate  cost  functionals  and  their  gradients.  We  do  a  comparison  with  a  problem 
considered  in  Chapter  11  of  [11].  In  this  chapter,  Burkardt  presents  what  he  refers  to  as 
a  discretized  sensitivity  failure.  We  will  evaluate  his  results  and  show  that  the  process  of 
adaptive  mesh  refinement  is  key  to  obtaining  good  cost  function  and  gradient  evaluations 
in  order  to  prevent  inaccurate  results  from  an  optimization  code. 

We  look  at  the  numerical  experiment  carried  out  in  §11.2  of  [11].  We  evaluate  the  following 
cost  functional: 


Jh{u{x,  y ),  q)  =  |  Y^{uh{ 3,  Vi)  -  uTaTg{ 3,  yi))2.  (5.1) 

i—1 

Here,  P  is  the  number  of  matching  points  and  for  the  results  we  present  herein  P  was  fixed 
at  15.  The  matching  points,  (3,  y *),  i  =  1, 2,  -  •  • ,  P,  were  evenly  distributed  along  the  line 
x  =  3,  y  £  [0,3].  Also,  the  target  profile,  uTaTg(3,y),  was  obtained  by  computing  a  finite 
element  solution  with  the  shape  parameter,  qTarg  =  (0.375, 0.5, 0.375)r,  and  A  =  0.5,  and 
then  interpolating  to  obtain  the  values  for  uT( 3,  y{):  i  =  1, 2,  •  •  •  ,  P. 

The  gradient  of  the  cost  function  with  respect  to  q  is  expressed  as 

V,j*  =  X>‘(3,  yi)  -  uTar9(3,  Vi)  As,  Si).  (5.2) 

The  SEM  uses  the  discrete  approximation  of  the  continuous  SE  as  before  to  approximate 
the  value  of  the  cost  gradient.  We  will  denote  this  approximations  of  the  cost  gradient  by 
Vq  Jh,  thus  we  have 


p 

Vq  *  Vq  Jh  =  £(u*(3,  yi)  -  uTar9( 3,  yi))sh(3,  Vi).  (5.3) 

*=  1 

In  order  to  nearly  duplicate  the  numerical  experiment  carried  out  by  Burkardt,  we  have 
L  =  10,  Re  =  1  and  we  generated  the  matching  profile  from  the  initial  mesh  for  qTorA  We 
calculated  the  values  of  the  cost  functional  Jh  and  its  gradient  along  a  line  parameterized 
by  5,  connecting  q  =  (-0.117, 0.419, -0.149)T  at  S  =  0  and  q  =  (0.375, 0.5, 0.375)r  at 
S  =  25.  This  is  the  same  line  along  which  Burkardt  explored.  His  results  are  presented  in 
Table  (11.2)  on  page  153,  and  Figures  (11.3)  and  (11.4),  page  145  of  [11]. 
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93 

Jh-1 03 

(VqJ^-lO3 

Jh-1 03 

(VqJft-5)-103 

Initial  Mesh 

03  Mesh 

SI 

0.419 

-0.149 

16.9 

-0.99 

12.1 

-0.56 

1$ 

0.422 

-0.128 

16.3 

-1.03 

11.7 

-0.57 

H 

0.425 

-0.107 

15.5 

-1.00 

11.3 

-0.58 

3 

-0.057 

0.428 

-0.086 

14.6 

-0.94 

10.8 

-0.60 

4 

-0.038 

0.432 

-0.065 

13.7 

-0.91 

10.3 

-0.61 

5 

-0.018 

0.435 

-0.044 

12.9 

-0.92 

9.95 

-0.62 

6 

0.001 

0.438 

-0.023 

12.2 

-1.04 

9.49 

-0.63 

7 

0.020 

0.441 

-0.002 

11.4 

-1.00 

9.02 

-0.65 

8 

0.040 

0.444 

0.018 

10.6 

-0.96 

8.56 

-0.67 

9 

0.060 

0.448 

0.039 

9.72 

-0.99 

8.00 

-0.68 

10 

0.079 

0.451 

0.060 

9.02 

-1.15 

7.48 

-0.70 

11 

0.099 

0.454 

0.081 

8.35 

-1.07 

6.95 

-0.72 

12 

0.119 

0.457 

0.102 

7.46 

-1.08 

6.39 

-0.74 

13 

0.138 

0.461 

0.123 

6.21 

-0.77 

5.76 

-0.75 

14 

0.158 

0.464 

0.144 

5.62 

-0.78 

5.17 

-0.76 

15 

0.178 

0.467 

0.165 

5.00 

-0.79 

4.56 

-0.77 

16 

0.197 

0.470 

0.186 

4.38 

-0.78 

3.94  1 

-0.76 

17 

0.217 

0.474 

0.207 

3.70 

-0.76 

3.27 

-0.75 

18 

0.237 

0.477 

0.228 

3.08 

-0.74 

2.66 

-0.72 

19 

0.256 

0.480 

0.249 

2.47 

-0.70 

2.07 

-0.68 

20 

0.276 

0.483 

0.270 

1.89 

-0.64 

1.51 

-0.62 

21 

0.296 

0.487 

0.291 

1.33 

-0.56 

0.99 

-0.53 

22 

0.316 

0.490 

0.312 

0.72 

-0.50 

0.57 

-0.42 

23 

0.335 

0.493 

0.333 

0.32 

-0.35 

0.25 

-0.29 

24 

0.355 

0.496 

0.354 

0.11 

-0.21 

0.05 

-0.13 

25 

0.375 

0.500 

0.375 

0.00 

0.00 

0.01 

0.05 

26 

0.394 

0.503 

0.396 

0.06 

0.16 

0.16 

0.28 

27 

0.414 

0.506 

0.417 

0.31 

0.38 

0.53 

0.53 

28 

0.434 

0.509 

0.438 

0.81 

0.63 

1.16 

0.81 

29 

0.453 

0.513 

0.459 

1.63 

0.92 

2.08 

1.13 

30 

0.473 

0.516 

0.480 

2.72 

1.22 

3.30 

1.48 

Table  5.3:  Values  of  Jh  along  a  Line,  Re  —  1,  and  A  =  0.5 


Dawn  L.  Stewart 


Chapter  5.  Numerical  Results  for  2-D  Problems 


112 


Our  results  are  displayed  in  Table  (5.3)  and  Figures  (5.30)  -  (  5.31).  Burkardt  reported 
a  local  maximum  at  S  =  5.  As  seen  in  Figure  (5.30),  we  do  not  get  the  same  results  on 
the  initial  mesh.  In  order  to  determine  if  the  difference  was  due  to  not  having  sufficient 
accuracy  on  the  initial  mesh,  we  also  presented  the  results  for  meshes  which  were  refined  on 
the  flow  field.  It  is  possible  that  the  difference  between  our  results  and  those  of  Burkardt 
could  be  due  to  our  having  fixed  the  value  of  A  at  0.5.  It  is  also  possible  that  the  difference 
is  a  result  of  numerical  inaccuracies  due  to  discretization  differences. 

We  did  one  more  numerical  experiment  exactly  similar  to  the  one  described  above  with 
one  exception,  we  now  used  Re  =  100.  Figure  (5.32)  -  (5.33)  show  the  cost  function  and 
gradient  along  the  same  line  explored  above.  Note  that  adapting  the  mesh  is  even  critical 
to  obtaining  accurate,  smooth  cost  and  gradient  approximations  at  this  Reynolds  number. 


Chapter  6 
Conclusions 


The  primary  goal  of  this  work  was  to  develop  and  analyze  new  algorithms  for  accurate 
computation  of  design  sensitivities.  We  focused  on  finite  element  schemes  applied  to  the 
continuous  sensitivity  equation.  The  SEM  is  more  efficient  than  standard  finite  difference 
schemes.  However,  accuracy  has  been  an  issue  for  flow  problems.  Our  work  concentrated 
on  two  methods  for  improving  accuracy  without  compromising  speed. 

The  first  method  used  enhanced  spatial  derivative  approximations  in  the  sensitivity  equa¬ 
tion.  Local  and  global  projection  techniques  were  employed  to  obtain  higher  accuracy  in 
the  derivative  approximations.  These  projection  techniques  not  only  provided  better  sen¬ 
sitivity  approximations  at  a  fixed  parameter  value,  but  also  stabilized  the  numerics  over 
greater  parameter  ranges.  For  the  1-D  model  problem,  the  local  projection  did  better  than 
the  global  projection  at  most  parameter  values  when  the  mesh  was  coarse.  As  the  mesh 
refined,  however,  the  use  of  global  projection  techniques  provided  more  accurate  sensitiv¬ 
ities.  At  this  point  there  are  no  theoretical  results  to  explain  this  behavior.  This  issue  is 
something  that  needs  to  be  addressed  in  future  work. 

The  second  technique,  mesh  refinement  based  on  errors  in  the  flow,  also  dramatically  im¬ 
proved  the  sensitivity  approximations.  This  was  especially  true  in  the  case  of  high  Reynolds 
number  calculations  or  in  shape  sensitivity  problems,  as  in  the  case  of  flow  over  a  bump.  In 
addition,  we  showed  that  at  higher  Reynolds  numbers  adapting  on  errors  in  the  sensitivities, 
as  well  as  the  flow,  improved  the  accuracy  of  the  sensitivity  approximations.  The  refine¬ 
ment  is  important  in  the  case  of  shape  problems  since  accurate  gradient  approximations  of 
the  state  are  needed  to  specify  the  sensitivity  boundary  conditions. 

We  discovered  that  when  these  two  techniques  were  used  for  approximating  sensitivities  for 
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use  in  the  calculation  of  cost  function  gradients,  the  convergence  properties  of  optimization 
algorithms  was  greatly  improved. 

There  are  many  issues  which  were  not  addressed  in  this  dissertation  which  will  be  the 
subject  of  future  research.  These  include  the  development  of  schemes  specifically  targeted 
to  improve  the  accuracy  of  local  projections  at  the  boundaries.  Also,  a  study  of  diagonal- 
ization  methods  for  use  in  optimization  algorithms  must  be  carried  out  to  determine  how 
the  trade-off  between  the  accuracy  of  cost  function/gradient  approximations  and  speed  of 
convergence.  Finally,  improved  error  estimates  are  needed  to  prove  convergence  of  these 
numerical  techniques  and  to  evaluate  rates  of  convergence  for  the  different  state  derivative 
approximations. 
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